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I.   INTRODUCTION 
A.      HISTORY 

Satellite  or  spacecraft  reentry  related  work  has  been  done  in  this  country  for  nearly 
five  decades.  With  the  advent  of  the  rocket  in  World  War  II,  the  military  saw  the  first 
application  of  spacecraft  reentry  prediction  as  the  ballistic  trajectory  calculation  of  the 
target  aimpoint.  Prior  to  this,  several  nations  had  "rocket  societies,"  however,  none  had 
actually  applied  the  use  of  "space"  to  a  technology  or  an  industry.  Thus,  the  first 
application  of  space  was  by  the  military  for  the  purpose  of  war.  [Ref.  l:p.  13] 

On  October  4,  1957,  the  former  Soviet  Union  changed  the  perception  of  space  held 

by  most  Americans  and  possibly  the  world.   It  was  the  launch  of  Sputnik  I  that  caused, 

then  Senator,  Lyndon  B.  Johnson  to  remark: 

That  sky  had  always  been  so  friendly,  and  had  brought  us  beautiful  stars  and 
moonlight  and  comfort;  all  at  once  it  seemed  to  have  some  question  marks  all  over 
it.  [Ref.  2:p.  13] 

When  the  former  Soviets  again  launched  another  rocket  barely  one  month  later  and 
sent  a  dog  into  space,  onboard  Sputnik  II,  the  world  now  saw  the  first  space  traveler  and 
the  dream  of  humans  in  space  became  more  real  than  dream  [Ref.  l:p.  20]. 

After  a  dismal  failure  of  the  United  States'  first  attempt  to  launch  a  satellite  into 
orbit  with  Vanguard  I,  which  was  dubbed  by  some  reporters  as  "Kaputnik"  after  it 
exploded  on  the  launch  pad  on  December  16,  1957,  Explorer  I  was  successfully  launched 


on  January  31,  1958.    The  space  race  had  begun,  and  both  countries  were  pushing  to 
have  the  first  human  in  space.  [Ref.  l:p.  20],  [Ref.  2:p.  16] 

With  the  goal  of  putting  humans  into  space  came  the  necessary  requirement  of 
providing  for  the  safe  return  of  those  humans  back  to  Earth,  unlike  the  first  dog  in  space 
which  perished  after  about  a  week  when  the  oxygen  supply  was  exhausted  [Ref.  2:p.  16]. 
This  was  the  beginning  of  the  most  intense  research  period  into  the  reentry  process  over 
the  entire  history  of  space  exploration. 

By  April  1972,  there  were  at  least  44  reported  instances  where  man-made  space 
objects  had  impacted  on  the  Earth  [Ref.  3:p.  383].  By  March  1978,  the  total  count  of 
man-made  objects  placed  into  Earth  orbit  was  10,791  [Ref.  4:p.  107].  By  Aug  1991,  the 
cumulative  count  of  objects  ever  placed  into  space  was  21,231  with  14,417 
decays/reentries,  leaving  6,814  objects  in  Earth  orbit  [Ref.  5:p.  v]. 

The  return  or  reentry  of  manned  spaceflights  has  been  covered  by  the  various 

media  sources  with  varying  degrees  of  intensity  based  upon  the  "newsworthiness"  of  the 

event  and  sensed  public  interest: 

Solar  Max  satellite  plunges  to  Earth.    [Ref.  6:p.  21(N)] 

Spacecraft's  study  of  sun  ends  today:    Solar  Max  heads  for  fiery  reentry.    [Ref. 
7:p.  A3] 

If  you  see  a  shooting  star  Dec.  8,  make  a  fast  wish  for  a  deep  cave.  [Ref.  8:p.  Bl] 

1.       Case  Studies 

It  is  useful  for  the  reader  to  understand  the  motivation  for  studying  the 
reentry  process.  Obviously,  the  routine  recovery  of  numerous  U.S.  and  Soviet  spacecraft 


must  imply  that  the  reentry  process  is  adequately  understood.  In  the  case  of  controlled 
reentry  and  space  vehicles  designed  for  reentry,  this  may  well  be  the  case;  however,  in 
the  uncontrolled  reentry  case,  this  is  not  as  easy  to  conclude. 

a.    Sputnik  TV 

September  5,  1962,  the  first  recorded  impact  of  a  man-made  space  object 
in  the  United  States  was  satellite  number,  1960  el,  Sputnik  IV.  This  satellite  reentered 
the  Earth's  atmosphere  over  North  America  and  many  fragments  of  the  reentering  debris 
ended  their  orbit  over  the  state  of  Wisconsin.  The  largest  fragment  of  which  was  found 
to  weigh  approximately  21  pounds  (9.49  kg)  impacted  in  the  city  of  Manitowoc, 
Wisconsin  at  approximately  0530  local  time.  [Ref.  9:p.  1] 

Sputnik  IV,  which  had  been  launched  on  May  14,  1960,  was  designed 
to  test  life-support  systems  for  Soviet  manned  space  flight.  On  May  19,  a  planned 
deorbit  maneuver  failed  and  the  pressure  vessel  separated  from  the  cabin,  leaving  the 
cabin  in  a  "lopsided"  orbit.  [Ref.  9:p.  2]  There  were  a  total  of  nine  separate  orbiting 
objects  associated  with  satellite  1960  el .  The  first  of  these  objects  to  reenter  was  the  last 
stage  rocket  body,  which  reentered  on  July  17,  1960.  By  July  1,  1961,  six  of  the  nine 
pieces  had  reentered.  The  payload  reentered  on  September  5,  1962  leaving  only  one 
other  object  in  orbit.  [Ref.  9:p.  3] 

The  U.S.  Space  Detection  and  Tracking  System  (SPADATS)  predicted 
Sputnik  IV  would  reenter  on  or  about  September  6,  1962.  Moonwatch  observers,  teams 
of  volunteers  who  attempted  to  observe  reentry  events  with  either  the  naked  eye  or 
telescopes,  were  notified  of  the  SPADATS  prediction  on  August  29,  1962.  [Ref.  9:p.  10] 


It  was  based  on  these  observations  that  the  fragments  of  Sputnik  IV  were  eventually 
located,  with  the  exception  of  the  largest  fragment,  which  was  found  by  a  routine  police 
patrol  of  the  city  of  Manitowoc.  Figure  1  shows  the  satellite  ground  trace  over  the  last 
one-half  revolution  [Ref.  9:p.  20].  Figure  2  shows  the  ground  trace  of  the  last  one 
hundred  nautical  miles  [Ref.  9:p.  12].  Figure  3  shows  the  impact  location  of  objects 
recovered  in  Manitowoc  [Ref.  9:p.  17].  The  objects  found  on  the  roof  of  the  church 
annex,  noted  in  Figure  3,  consisted  of  15  small  spherules  approximately  1/8  inch  in 
diameter.  These  were  about  325  feet  further  downrange  from  the  initial  mass  which 
impacted  near  the  intersection  of  Park  St.  and  8th  St. 

There  was  no  loss  of  life  and  the  only  reported  property  damage  was  the 
impact  impression  left  in  the  street. 

b.    Skylab 

In  May  1973,  the  last  U.S.  Saturn  V  rocket  launched  Skylab  into  orbit 
approximately  270  nautical  miles  (nm)  above  the  Earth.  Skylab  was  the  first  U.S.  space 
station  and  the  centerpiece  of  the  U.S.  space  program  since  the  last  U.S.  moon  mission. 
Skylab  would  remain  in  orbit  until  July  11,  1979.  [Ref.  10:p.  1] 

Knowing  the  inevitability  of  decay  and  reentry  of  low  Earth  orbits, 
NASA  contracted  the  Lockheed  Missile  and  Space  Company  (LMSC)  to  investigate  the 
predicted  reentry  and  breakup  of  Skylab  three  years  prior  to  its  launch.  This  initial 
report  predicted  that  Skylab  would  begin  to  breakup  at  an  altitude  of  65  nm  (120  km)  and 
that  debris  would  fall  3600  nm  downrange  from  the  initial  breakup.  [Ref.  10:p.  2] 
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Figure  1:   Sputnik  IV  Final  Revolution 
[Ref .  9] 


Figure  2:   Sputnik  IV  Impact  Area/Debris  Location 
[Ref.  9] 
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Figure  3:    Sputnik  IV  Impact  Area  Ground  Trace 
[Ref.  9] 


Based  on  the  post-flight  reconstruction  of  the  reentry  of  Sky  lab,  Marshall 
Space  Flight  Center  (MFSC)  was  able  to  conclude  that  telemetry  signals  were  still  being 
sent  from  Skylab  as  it  passed  over  Bermuda  and  Ascension  Islands.  Also,  each  of  these 
tracking  stations  reported  only  one  radar  contact.  It  was  therefore  concluded  that  Skylab 
was  intact  and  breakup  had  not  begun  although  the  altitude  had  decreased  to  57  nm.  [Ref. 
10:p.  5]  Survivability  underestimation  has  been  typical  of  past  analysis  as  will  be 
discussed  in  Chapter  III. 

Figure  4  shows  the  ground  trace  of  Skylab  over  the  last  one-quarter 
revolution  and  denotes  the  impact  footprint  (debris  dispersion  area)  [Ref.  10:p.  18]. 
Table  1  is  a  partial  listing  of  recovered  debris  from  Skylab  and  the  recovery  location 
[Ref.  10:p.  10]. 

Again,  there  was  no  loss  of  life,  no  human  injury  and  no  property 
damage  as  a  result  of  this  uncontrolled  reentry. 

c.     Cosmos  954 

January  24,  1978  marked  the  first  uncontrolled  reentry  and  Earth  impact 
of  a  nuclear  powered  artificial  satellite  [Ref.  3:p.  384].  This  was  also  the  first  case  of 
"significant"  property  damage  due  to  an  artificial  Earth  satellite  impact.  It  was  reported 
that  Canada  spent  over  1 1  million  dollars  and  the  United  States  spent  nearly  3  million 
dollars  in  the  location  and  recovery  of  radioactive  debris  [Ref.  3:p.  386].  Figure  5 
shows  the  reentry  ground  trace  and  impact  dispersion  area  of  Cosmos  954  [Ref.  ll:p. 
303]. 


Figure  4:    Cosmos  954  Impact  Area  /  Debris  Dispersion 
[Ref.  10] 


Table  I: 
[Ref.  10] 


SKYLAB  RECOVERED  DEBRIS  (PARTIAL  LISTING) 


ITEMS 

PROBABLE  SOURCE 

LOCATION 

Charred  Fragments 

OWS 

33. 9S, 

121. 9E 

(In  Esperance) 

Burned  Material 

ows 

33. 9S, 

121. 9E 

(In  Esperance) 

Aluminum  356  Casting 

OWS 

33. 7S, 

122. IE 

(20  mi  NE  of 
Esperance) 

Foam  Fiberglass 
Beam  Section 

OWS 

33. 9S, 

122. 0E 

(9  mi  E  of 
Esperance) 

H20  Tank 
Aft  End 

OWS 

33.8S, 

122. 0E 

(9  mi  NE  of 
Esperance) 

H20  Tank 

OWS 

33. 9S, 

122. IE 

(10  mi  E.  of 
Esperance) 

10'  Steel  Strip  ' 

OWS 
(H20  Tank) 

33. 9S, 

122. 3E 

(25  mi  E  of 
Esperance) 

Heat  Exchanger 

OWS 
(H20  Cooler) 

33. 9S, 

122. IE 

(12  mi  E.  of 
Esperance) 

Segment  of 
Fiberglass  Sphere 

OWS 

33. 9S, 

122. IE 

(11  mi  E  of 
Esperance) 

Insulation 

OWS 
(Bulkhead) 

33. 9S, 

122. IE 

(11  mi  E  of 
Esperance) 

Aluminum  Gear 
and  Housing 

OWS 
(Urine  Separator) 

33. IS, 

122. 5E 

(40  mi  NE  of 
Esperance) 

N2Tank 

AM 

33. 2S, 

122. 6E 

(60  mi  NE  of 
Esperance) 

Electronics  Module 

AM 

33. 5S, 

122.3E 

(35  mi  NE  of 
Esperance) 

N-Sphere 

AM 

33. 5S, 

122. 8E 

(49  mi  ENE  of 
Esperance  in 
Neridup  area) 

Pressure  Tank 

IU 

33. 2S, 

122. 6E 

(60  mi  NE  of 

Esperance) 
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Figure  5:   Skylab  Final  One-Quarter  Revolution 
[Ref.  11] 
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The  Outer  Space  Treaty  of  1967,  the  Liability  Convention  of  1972  and 
the  Registration  Convention  of  1976  are  the  three  principal  documents,  drafted  by  the 
United  Nations  Committee  on  the  Peaceful  Uses  of  Outer  Space  and  its  Legal 
Subcommittee,  which  define  ownership  responsibilities  of  space  objects  [Ref.  12:p.  457]. 
Although  the  Outer  Space  Treaty,  Article  VIII,  provides  in  explicit  terms  the 
requirements  that  the  registry  state  (country  of  origin)  retains  jurisdiction  and  control 
over  its  satellites  while  in  space  or  on  a  celestial  body,  [Ref.  3:p.  389]  there  may  be  no 
legal  requirement  under  the  Liability  Convention  or  the  Rescue  and  Return  Agreement 
for  the  former  Soviet  Union  to  reimburse  either  the  U.S.  or  Canada  for  the  cleanup  cost 
of  Cosmos  954  [Ref.  3:p.  387]. 

It  is  because  of  the  potential  for  the  loss  of  human  life,  significant 
property  damage  and  legal  responsibility  based  in  international  law  that  the  reentry  of 
uncontrolled  artificial  satellites  must  be  further  investigated.  In  the  case  of  Cosmos  954, 
on  the  morning  of  reentry,  the  Soviets  predicted  reentry  impact  near  the  Aleutian  Islands 
(52°N,173°W)  and  the  North  American  Air  Defense  Command  (NORAD)  predicted 
reentry  impact  near  Hawaii  (19°N,156°W)  [Ref.  4:p.  110].  The  physics  of  uncontrolled 
reentry  and  the  modeling  of  the  reentry  process  must  be  better  understood  in  order  to 
improve  the  accuracy  of  reentry  predictions. 

B.      PROBLEM  STATEMENT 

The  probability  of  space  objects  reentering  the  Earth's  atmosphere  and  surviving 
to  impact  is  relatively  small;  most  artificial  objects  burn  up  in  the  atmosphere  before  they 
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impact  [Ref.  13:p.  1].  Falling  debris  nonetheless  presents  a  potential  hazard  to  people 
and  property.  The  most  difficult  problem  associated  with  reentry  impact  prediction  is 
not  to  determine  what  will  survive  to  impact  in  as  much  as  where  impact  will  occur. 
[Ref.  5:p.  v] 

The  significant  factors  related  to  the  prediction  of  time  to  impact  and  location  of 
the  reentry  are  a  lack  of  observational  data  over  the  entire  orbit  and  a  lack  of  "precise" 
mathematical  models  which  accurately  describe  the  physical  processes  occurring  during 
reentry.  The  major  parameters  which  contribute  to  the  uncertainties  of  the  reentry 
prediction  are:  [Ref.  5:p.  v] 


1.  Atmospheric  density  variations—atmospheric  density  is  strongly  influenced  by 
solar  and  geomagnetic  activity,  both  of  which  are  difficult  to  forecast. 

2.  Aerodynamic  force  models-aerodynamic  forces  are  a  function  of  attitude,  lift 
and  drag  coefficients,  gas-surface  interactions  and  gas  dynamics  such  as 
continuum  flow  or  free-flow  regimes. 

3.  Spacecraft  attitude  motion~the  attitude  of  the  object  and  how  it  changes  with 
time  is  an  important  factor  in  estimating  the  aerodynamic  forces  encountered 
during  reentry. 

4.  Changes  in  configuration— ablation  and  fragmentation  cause  changes  in 
configuration  (profile  area  and  mass/area  loss)  which  can  change  the 
aerodynamic  forces  experienced  either  as  an  increase  or  decrease  in  net  force. 


The  lack  of  regularly  spaced  observational  data  over  the  entire  orbit  can  severely 
handicap  the  efforts  to  predict  reentry,  especially  in  the  final  phase.  Since  most 
observational  data  is  from  radar  tracking  stations  and  the  object  is  in  a  very  low  orbit, 
the  time  in  view  is  of  short  duration  as  well  as  limited  by  the  geographic  location  of  the 
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tracking  stations.  The  availability  of  good  quality  tracking  data  is  required  for  accurate 
reentry  predictions  and  could  compensate  for  inherent  deficiencies  in  the  models.  [Ref. 
5:p.  vi] 

C.      PURPOSE 

This  research  was  initiated  by  the  Air  Force  Space  Command.  The  purpose  of  this 
thesis  is  threefold: 


1.  Conduct  a  comprehensive  literature  survey  in  the  area  of  artificial  satellite 
reentry  specifically,  uncontrolled  orbit  decay  and  reentry  into  the  Earth's 
atmosphere. 

2.  Describe  the  "state-of-the-art"  of  reentry/impact  prediction  techniques. 

3.  Define  critical  areas  of  research  where  increased  emphasis  is  required  in  order 
to  improve  the  accuracy  of  the  reentry  prediction. 


The  focus  of  this  thesis  is  the  physical  processes  of  the  uncontrolled  reentry  from 
the  final  few  revolutions  to  impact.  It  is  not  the  intended  purpose  of  this  thesis  to 
examine  in  detail  or  focus  on  the  following  aspects  of  the  reentry  prediction  problem: 

1.  Atmospheric  density   models 

a.  Solar  activity  and  influence  on  reentry 

b.  Geomagnetic  influence  on  reentry 

2.  Long-term  orbital  lifetime  prediction 
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It  is  necessary,  however,  to  discuss  these  aspects  of  the  reentry  process  in  order  to 
completely  study  the  physical  phenomena  associated  with  reentry.  These  topics  are  the 
subject  matter  of  Chapter  II  and  are  considered  the  fundamental  background  material 
required  for  a  further,  more  detailed  study  of  the  reentry  process. 

D.      METHODOLOGY 

The  primary  research  methodology  of  this  thesis  was  a  comprehensive  literature 
survey  of  military  and  civil  aerospace  data  bases.  The  literature  search  of  Department 
of  Defense  (DoD)  data  bases  was  restricted  to  unclassified  work.  The  literature  search 
was  conducted  at  or  through  the  following  activities: 

1.  Naval  Postgraduate  School,  Dudley  Knox  Library 

2.  Hanscom  Air  Force  Base  Research  Library 

3.  University  of  Colorado  at  Boulder,  Astrophysics  and  Engineering  libraries 

4.  AFSPACECOM  Astrodynamics  Division  (CNY)  technical  library 

5.  United  States  Air  Force  Academy  Library 


The  secondary  research  methodology  of  this  thesis  was  personal  interviews  with 
"experts"  in  the  study  of  reentry  or  reentry  related  fields.  The  activities  which  were 
contacted  or  visited  personally  include: 

1.  TRW  Corporation 

2.  The  Aerospace  Corporation 
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3.  Phillips  Laboratories 

4.  Air  Force  Geophysics  Laboratory  (AFGL) 

5.  Naval  Space  Surveillance  Center  (NAVSPASUR) 

6.  Air  Force  Space  Command  (AFSPACECOM) 

7.  United  States  Space  Command  (USSPACECOM) 


The  primary  data  bases  and  periods  over  which  the  literature  survey  was  conducted 
include: 

1.  NASA  1940-1993 

2.  DTIC  1950  -  1993 

3.  IAA  1960  -  1993 

4.  STAR  1960  -  1993 

5.  DIALOG  1970  -  1993 


These  data  bases  were  searched  with  similar  strategies  in  an  effort  to  pull  as  many 
"original"  articles  as  possible,  while  concurrently  verifying  the  search  process  by 
producing  numerous  "duplicate"  articles  found  in  other  data  bases.  The  primary  search 
terms  included:  reentry,  reentry  prediction,  satellite  reentry,  atmosphere  reentry, 
spacecraft,  spacecraft  reentry,  uncontrolled  reentry,  reentry  impact  prediction,  satellite, 
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descent  trajectory,  atmosphere  drag,  atmospheric  density  models,  satellite  orbit  decay, 
reentry  heating,  reentry  dynamics  and  reentry  effects. 

In  order  to  prevent  confusion  of  common  variables  when  derived  in  multiple 
sources,  the  equations  and  variable  nomenclature  presented  throughout  this  thesis  are  as 
presented  in  the  original  works,  with  the  exception  of  minor  changes  made  as  noted.  An 
example  of  such  a  change  is  the  flight  path  angle,  7,  presented  in  Chapter  II.  Since  this 
is  a  survey  of  the  literature,  which  dates  back  to  the  1950's,  the  authors  have  made  a 
conscious  effort  to  preserve  the  flavor  of  the  individual  works  and  no  attempt  has  been 
made  to  standardize  the  nomenclature.  However,  this  is  an  object  for  further  study  as 
indicated  in  Chapter  VI. 
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H.      FUNDAMENTALS  OF  ATMOSPHERIC  REENTRY 

A.      TYPES  OF  REENTRY 

The  reentry  trajectories  of  space  vehicles  can  be  classified  into  two  types,  either 
uncontrolled  or  controlled.  The  focus  of  this  thesis  is  to  investigate  the  uncontrolled 
reentry  of  satellites,  however,  three  major  types  of  controlled  reentry,  ballistic,  gliding 
and  skip,  as  shown  in  Figure  6,  will  be  described  for  completeness  [Ref.  14:p.  6]. 

1.  Uncontrolled  Reentry 

Uncontrolled  reentry  may  be  the  result  of  an  unrecoverable  satellite  subsystem 
failure  or  the  end  of  the  satellite's  operational  life.  The  flight  path  angle  is  usually  much 
less  than  one  degree  and  lift  is  considered  negligible.  Typically,  very  large  uncertainties 
in  impact  point  predictions  are  created  by  decay-induced  uncontrolled  reentries.  [Ref. 
15:p.  44] 

2.  Controlled  Reentry 

During  a  controlled  reentry,  the  vehicle's  aerodynamic  and  heating  loads  are 
maintained  within  acceptable  limits  by  controlling  the  effects  of  lift  and  drag  forces  on 
the  vehicle  throughout  the  flight.  This  is  accomplished  through  a  carefully  designed 
space  vehicle,  flight  trajectory  and  possibly  a  precision  guidance  system.  Controlled 
reentry  spans  an  aerodynamic  flight  regime  from  subsonic  to  Mach  25  and  beyond  for 
hyperbolic  reentry.  [Ref.  16:p.  231] 
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Decaying  orbit 
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Figure  6:   Types  Of  Reentry 
[Ref.  14] 
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a.  Ballistic  Reentry 

A  ballistic  reentry  is  characterized  by  sufficiently  steep  reentry  angles 
where  the  force  of  lift  is  assumed  to  be  negligible  [Ref.  14:p.  2].  The  ability  to  control 
the  reentry  velocity,  flight  path  angle,  ballistic  coefficient  and  atmospheric  properties 
determines  the  accuracy  of  ballistic  reentry  vehicles.  Impact  accuracies  within  the 
intended  target  vary  from  20  km  for  shallow  angle  reentries,  such  as  Mercury  type 
vehicles,  to  an  accuracy  of  a  few  hundred  meters  for  a  steep  angle  reentry  of  an 
intercontinental  ballistic  missile  type  vehicle.  [Ref.  16:pp.  237-242] 

b.  Gliding  Reentry 

A  gliding  reentry  is  characterized  by  a  glide  slope  rather  than  a  reentry 
trajectory.  During  a  gliding  reentry,  a  vehicle  such  as  the  space  shuttle  creates  enough 
lift  to  maintain  a  long  hypersonic  glide  at  a  small  flight  path  angle  [Ref.  16:p.  242].  A 
measure  of  the  vehicle's  lift  that  influences  the  descent  path  and  cross  range  capability 
of  a  vehicle  is  the  lift-to-drag  (L/D)  ratio  [Ref.  15:p.  47].  By  adjusting  the  vehicle  L/D 
ratio  or  bank  angle,  a,  (the  angle  between  the  vehicle  lift  vector  and  the  plane  containing 
the  vehicle  position  and  velocity  vectors),  the  extent  of  the  range  and  cross-range  can  be 
controlled  [Ref.  15:p.  44]. 

c.  Skip  Reentry 

A  skip  reentry  is  characterized  by  a  vehicle  whose  L/D  is  greater  than 
zero.  Sufficient  lift  is  produced  to  dominate  the  gravitational  and  centrifugal  forces.  If 
this  lift  is  combined  with  a  large  enough  initial  angle  of  descent,  a  reentry  trajectory  with 
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one  or  more  skips  may  be  produced.  As  the  vehicle  begins  to  reenter,  it  reaches  a 
minimum  altitude  where  it  begins  to  "pull-up"  due  to  the  lift  force  dominating  over 
gravity.  Eventually,  the  vehicle  will  exit  the  atmosphere  at  a  reduced  velocity.  If  the 
exit  velocity  and  flight  path  angle  are  correctly  controlled,  the  vehicle  will  achieve  a  brief 
orbital  phase  followed  by  a  second  reentry  downrange  from  the  first.  [Ref.  16:pp.  245- 
246] 

B.      MAJOR  FACTORS  INFLUENCING  UNCONTROLLED  REENTRY 

It  is  the  focus  of  this  thesis  to  describe  the  physical  processes  of  the  final  few 
revolutions  of  a  reentry  body's  orbit  and  the  reentry  phase  to  impact.  Reentry  trajectory 
techniques  differ  from  orbit  determination  techniques.  The  reentry  phase  is  very 
dynamic  as  opposed  to  the  exoatmospheric  phase  of  orbital  motion.  Specifically,  the 
reentry  phase  is  characterized  by: 

1.  Rapidly  decreasing  altitude 

2.  Rapidly  increasing  aerodynamic  heating  effects 

3.  Rapidly  changing  aerodynamic  load  effects 


The  reentry  phase  can  be  characterized  by  that  portion  of  the  trajectory  prior  to 
breakup  and  after  breakup.  Prior  to  breakup,  the  reentry  trajectory  equations  of  motion 
include  parameters  such  as  gravity,  atmospheric  density,  ballistic  coefficient,  position  and 
velocity.     Solutions  to  the  reentry  equations  of  motion  can  be  found  using  analytical, 
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semi-analytical  and  numerical  techniques.  Breakup  can  be  described  as  that  point  in  the 
trajectory  where  heating  effects  and  load  effects  cause  the  reentry  body  to  lose  its 
structural  integrity. 

1.       Modeling  Reentry  Equations  Of  Motion 

a.    Basic  Equations  of  a  Rigid  Body 

The  three  basic  equations  of  a  rigid  body  are 


*  =  V  (D 

dt 


?!£  =  F  (2) 

dt 


¥l=M  (3) 

dt 

where 

H       =  reentry  body's  angular  momentum  vector  relative  to  its  center  of  mass 

K       =  reentry  body's  linear  momentum  vector 

M      =  total  moment  force  vector  relative  to  the  center  of  mass 

F        =  total  force  vector  acting  on  the  body 

r        =  position  vector  of  the  body 

V       =  velocity  vector  of  the  center  of  mass 

Equations  (1)  and  (2)  are  the  kinematic  and  force  equations,  respectively. 
When  these  equations  are  coupled,  they  yield  one  second-order,  non-linear,  vector 
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differential  equation  which  describes  three-dimensional  motion  of  the  reentry  body's 
center  of  mass,  given  by 

—  =  g+a,  (4) 

where  g  is  the  gravitational  field  and  adn%  represents  the  atmospheric  drag  acceleration. 
It  is  atmospheric  drag  that  causes  an  artificial  satellite  to  decay  and  reenter. 

Equation  (3)  describes  the  motion  of  a  body  about  its  center  of  mass, 
commonly  referred  to  as  attitude  motion.  The  reentry  body's  attitude  is  related  to  angle 
of  attack,  a,  and  is  an  important  parameter  which  will  be  discussed  in  the  next  section. 

When  the  system  of  equations  (1)  through  (3)  are  coupled  along  with  the 
attitude  parameters,  it  yields  two  second-order,  non-linear,  vector  differential  equations 
or  six  second-order  scalar  differential  equations.  These  six  differential  equations 
completely  describe  the  six-dimensional  motion  of  and  about  the  center  of  mass.  This 
system  of  equations  is  commonly  referred  to  as  the  six-degree-of-freedom  model  which 
will  be  described  further  in  Chapter  III. 

b.    Modeling  Gravity 

In  equation  (4)  the  gravitational  field,  g,  may  be  written  as  an  inverse 
square  relationship 


g.-JLr  (5) 

r3 
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where 

\k  =  gravitational  constant  (Earth  =  3.9865  x  105  kmVsec2) 
or  more  accurately  via  the  geopotential  model.  The  geopotential  model  divides  the  Earth 
into  three  sets  of  geographically  divided  regions  described  by  latitude  and  longitude  as 
shown  in  Figure  7  [Ref.  17:p.  233].  The  significance  of  the  geopotential  model  is  its 
ability  to  describe  the  gravitational  field  more  accurately  than  a  point  mass  model  for  the 
Earth.  For  example,  in  the  region  of  reentry  below  120  km,  the  first-order  zonal 
harmonic  J2,  may  attain  a  magnitude  approaching  that  of  atmospheric  drag  [Ref.  18:p. 
40]. 

c.     Modeling  Atmospheric  Drag 

As  previously  mentioned,  the  force  that  causes  an  artificial  satellite's 
orbit  to  decay  is  atmospheric  drag.  The  atmospheric  drag  acceleration  vector,  in 
equation  (4),  acts  in  the  direction  opposite  of  the  satellite's  relative  velocity  vector  and 
is  given  by  [Ref.  17:p.  258] 


.-,  ■  -^p* 


where 

CD  =  satellite  drag  coefficient 

A  =  aerodynamic  effective  cross-section  area 

m  =  satellite  mass 

p  =  local  neutral  density  of  the  atmosphere 
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Figure  7:    Earth  Geopotential  Model 
[Ref.  17] 
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VA      =  I  V  -  we  x  r  I    (satellite's  airspeed) 
v        =  unit  vector  in  direction  of  (V  -  we  x  r) 
we       =  2t  rad/day 

V       =  satellite's  inertial  velocity  vector 
The  ballistic  coefficient,  from  equation  (6),  is  defined  as 

B«f2l  (7) 

m 

Atmospheric  drag  dissipates  the  satellite's  energy  which  in  turn  causes 
a  decrease  in  the  semi-major  axis.  Drag  affects  high  eccentricity  orbits  by  gradually 
decreasing  the  apogee  altitude,  while  maintaining  a  nearly  constant  perigee  altitude, 
resulting  in  circularizing  the  orbit  as  shown  in  Figure  8  [Ref.  17:p.  258].  This 
contraction  will  continue  until  the  satellite  begins  the  reentry  phase.  Under  certain 
simplifying  assumptions,  such  as  ignoring  the  rotation  of  the  atmosphere,  the  analytical 
results  describing  the  change  per  one  revolution  in  semi-major  axis  and  eccentricity  for 
orbital  decay  are  given  by  [Ref.  19:p.  230] 


Cyi^2r2^(l+gcos£)2 
(1  -ecosE) 


2 


Ae  =  -— 


m 

where 


r2n  (\+ecosE) 
Jo     \  1  -ecosE; 


'(cos  E+e)  dE 


(9) 


=  semi-major  axis 
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Figure  8:    Drag  Effect  On  High  Eccentricity  Orbits 
[Ref.  17] 
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e        =  eccentricity 

E        =  eccentric  anomaly 

(1)  Modeling  Atmospheric  Density.  The  Earth's  atmosphere  is 
primarily  composed  of  nitrogen  and  oxygen.  Solar  radiation  affects  the  dynamic 
properties  of  this  medium  by  constantly  changing  the  temperature,  pressure,  chemical 
constituents,  particulate  presence  and  electrical  properties  [Ref.  20:p.  5].  The  inability 
to  model  the  atmospheric  neutral  density  is  an  integral  part  of  the  satellite  reentry 
problem.  Neutral  density  is  defined  as  the  density  of  the  neutrally  charged  constituents 
of  the  atmosphere.  The  atmospheric  density  is  not  precisely  known  along  the  satellite 
path  because  it  varies  with  geographic  location,  solar  and  geomagnetic  conditions, 
altitude  and  time  [Ref.  17:p.  257]. 

Atmospheric  density  models  are  categorized  as  either  "theoretical," 
or  "empirical"  models.  Empirical  models  describe  the  phenomena  of  the  Earth's 
atmosphere,  based  on  a  summary  of  observed  data  and  are  constructed  independently  of 
the  laws  of  physics.  Theoretical  models  apply  the  laws  of  physics.  Empirical  and 
theoretical  models  may  have  overlapping  domains.  For  instance,  data  used  to  construct 
a  particular  empirical  model  may  be  smoothed  and/or  extrapolated  based  on  theoretical 
considerations.  Likewise,  theoretical  models  must  be  compared  to  empirical  data  in 
order  to  define  empirical  parameters  such  as  average  sea  level  values  of  pressure  and 
temperature  and  to  establish  assumption  limits.  [Ref.  21:pp.  B-l,B-2]  For  example,  the 
well  known  Jacchia  models  and  the  global  dynamic  models  are  semi-theoretical  and  semi- 
empirical  in  nature.    The  global  dynamic  models  describe  the  physical  and  chemical 
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processes  of  the  Earth's  coupled  thermosphere-ionosphere  system.  Solar  radiation  and 
auroral  input  measurements  are  used  in  the  physical  model  to  predict  the  time-dependent 
density  response.  [Ref.  22:p.  1] 

A    simplified   analytical    atmospheric    model   demonstrates    the 
relationships  between  temperature,  altitude  and  density 

P  -  p/*)  (10) 

where 

p0       =  sea  level  air  density 

z        =  altitude 

H       =  atmospheric  scale  height  (RT0  /  g) 

T0      =  atmospheric  temperature  at  sea  level 

R  =  gas  constant  (air) 
Equation  (10)  is  a  simple  analytic  relationship  where  the  atmospheric  density  decreases 
exponentially  with  altitude  when  gravity,  temperature  and  chemical  composition  are 
assumed  to  be  constant  at  all  altitudes.  This  model  represents  a  rough  approximation  of 
the  atmospheric  density.  Additional  physical  properties  and  levels  of  sophistication  can 
be  incorporated  into  this  simplified  model  in  order  to  better  describe  the  actual 
atmosphere. 

Empirical  atmospheric  density  models  use  a  variety  of  functions  to 
describe  the  atmosphere.    Some  of  these  functions  are  common  to  most  models.    The 
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following  list  provides  a  brief  description  of  these  common  functions:  [Ref.  21  :pp.  B-13-- 
17] 

1 .  Altitude—density  decreases  with  increasing  altitude. 

2.  Early  afternoon  bulge—maximum  daytime  temperatures  cause  the  density  to 
increase  at  satellite  altitudes. 

3.  Average  10.7  cm  (F,07)  flux— the  average  solar  power  per  unit  area  at  a 
frequency  of  2800  MHz  (X= 10.7  cm)  measured  at  various  Earth  locations.  The 
10.7  cm  flux  is  closely  correlated  to  the  extreme  ultra-violet  (EUV)  radiation 
which  heats  up  the  upper  atmosphere  and  is  used  as  an  indicator  of  the  EUV 
flux,  since  EUV  flux  is  absorbed  at  higher  altitudes  and  is  converted  to  heat. 

4.  Daily  10.7  cm  flux-accounts  for  the  rapidly  changing  values  of  the  F10.7 
measurement. 

5.  Geomagnetic  index  (A^-related  to  the  activity  of  charged  particles.  Usually 
modeled  as  a  correction  to  the  atmospheric  temperature. 

6.  Winds- speeds  up  to  300  m/s  have  been  observed  in  the  upper  atmosphere. 
Wind  can  significantly  change  the  drag  experienced  by  the  satellite  since  drag 
is  proportional  to  the  square  of  the  velocity  with  respect  to  the  surrounding  air. 


Numerous  empirical  atmospheric  density  models  have  been 
developed  since  the  launch  of  Sputnik  I  [Refs.  23-31].  Early  models  such  as  the  Jacchia 
70  and  71  were  derived  from  the  analysis  of  satellite  drag.  These  models  identified  the 
upper  atmosphere  as  dependent  on  solar  flux,  geomagnetic  index,  diurnal,  monthly  and 
seasonal  variations.  A  later  model,  the  Jacchia  77  model  incorporated  composition  data 
of  nitrogen  (Nj)  and  mono-atomic  oxygen  (O)  which  was  observed  from  satellite  mass 
spectrometers.    The  mass  spectrometer  and  incoherent  backscatter  (MSIS)  model  was 
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constructed    using    composition    data    derived    from    satellite    mass    spectrometers, 
accelerometers  and  ground  based  incoherent  backscatter  measurements. 

Several  comparison  studies  have  been  conducted  with  various 
empirical  atmospheric  density  models  to  determine  their  overall  accuracy  and  efficiency. 
These  studies  indicate  the  following  limitations  and  deficiencies  of  empirical  models 
[Refs.  21,  32-34] 


1.  Accuracies  of  models  have  remained  relatively  unchanged  for  the  past  two 
decades. 

2.  Statistical  analysis  of  measured  satellite  accelerometer  density  data  as  compared 
with  atmospheric  density  model  mean  values  and  standard  deviations,  indicate 
mean  value  accuracies  of  approximately  ±10%  with  standard  deviations  of 
approximately  ±15%. 

3.  Some   models  are  significantly  less  efficient  in  terms  of  computational  time. 

4.  The  F107  cm  measurement  does  not  adequately  represent  the  complex  interaction 
between  the  EUV  flux  and  the  thermosphere. 

5.  The  geomagnetic  index,  Ap,  or  3  hour  Kp  does  not  necessarily  represent  the 
physical  mechanism  that  causes  the  variation  in  atmospheric  density. 

6.  New  model  parameters  such  as  the  precipitation  index  (used  as  an  indicator  of 
magnetospheric  activity),  are  needed  to  model  the  real  physical  variations  in  the 
atmosphere. 


Additionally,  our  ability  to  predict  the  solar  flux  and  the 
geomagnetic  field  is  usually  difficult  and  unreliable  at  best.  As  a  consequence,  accurate 
forecasting  of  atmospheric  density  into  the  future  is  limited.      This  affects  the 


31 


determination  of  the  predicted  orbital  decay  rates  because  of  the  dynamic  dependence  on 
the  variations  of  the  F107  and  Ap  indices.  [Ref.  35 :p.  1] 

Considerable  progress  has  been  made  over  the  past  decade  in  the 
development  of  a  dynamic  atmospheric  model  of  the  coupled  thermosphere-ionosphere 
system.  The  National  Center  for  Atmospheric  Research  (NCAR)  has  developed  a  model 
as  the  thermosphere/ionosphere  general  circulation  model  (TIGCM).  The  pressure 
coordinate  primitive  equations  of  the  lower  atmospheric  meteorology  are  the  foundation 
of  the  TIGCM.  TIGCM  uses  20  minutes  of  CRAY-Y-MP  8/64  computer  time  per 
simulated  model  day  to  compute  the  prognostic  thermodynamic,  eastward  and  northward 
momentum  equations  and  diagnostic  equations  of  state  and  continuity.  The  model  utilizes 
a  5°  lat-long  grid  with  24  constant  pressure  surfaces  distributed  from  an  altitude  of  95 
to  500  km.  Density  data  collected  from  a  200  km  altitude  satellite  was  compared  with 
the  TIGCM  calculated  density  and  the  results  showed  that  the  TIGCM  calculations  were 
within  ±  9-12%  on  a  point  by  point  basis  along  the  satellite's  track.  A  new  version, 
TIE-GCM,  added  an  interactive  dynamo  model  to  calculate  electrodynamic  interactions 
between  the  thermosphere  and  the  ionosphere.  Results  indicate  this  version  model  is  able 
to  provide  improved  determination  of  thermospheric  density,  especially  during  disturbed 
geomagnetic  conditions.  [Ref.  22:pp.  1-6]  An  operational  version  of  this  model,  using 
a  vector  spherical  harmonic  (VSH)  technique  is  under  development  at  the  University  of 
Michigan.  With  further  development,  this  dynamic  atmospheric  model  may  be  able  to 
forecast  global  atmospheric  density  values  with  errors  less  than  10%.  [Ref.  36] 
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(2)  Modeling  Aerodynamic  Coefficients.  The  aerodynamic  coefficients, 
depend  upon  the  following  factors: 

1.  The  shape  and  dimensions  of  the  vehicle 

2.  The  vehicle  orientation  to  the  on-coming  air  flow  (angle  of  attack) 

3.  The  temperature  and  composition  of  the  neutral  atmosphere 

4.  The  gas-surface  interaction  phenomenon 

The  flow  phenomena  encountered  around  complex  shapes  further  complicates  and  varies 
this  parameter  in  the  aerodynamic  regimes:  free  molecular  flow,  transitional  flow,  and 
continuum  flow  [Ref.  37:p.  5].  A  CD  of  2.2  is  typically  used  as  a  constant  value  in  the 
free  molecular  flow  region  above  120  km.  However,  CDis  a  function  of  angle  of  attack, 
shape  and  flow  regime.  Therefore,  knowing  the  vehicle's  motion  about  its  center  of  mass 
is  a  critical  factor  in  determining  the  aerodynamic  forces  acting  on  the  vehicle.  These 
forces  dictate  the  orbit  decay  rate,  reentry  and  impact  location.  CD  can  vary  from  2  in 
the  free  molecular  flow  regime  to  a  value  much  less  than  1  in  the  continuum  flow 
regime. 

(3)  Ratified  Gas  Dynamics.  The  pressure  distribution  and  subsequent 
forces  imparted  by  the  near  flowfield  on  a  reentry  body  determines  the  forces  and 
moments  acting  on  the  body  via  the  aerodynamic  coefficients  [Ref.  20:p.  203].  There 
are  basically  five  flow  regimes  which  have  distinguishable  characteristics.  These  five 
flow  regimes  may  be  given  quantitative  definition  using  the  Knudsen  number  (Kn).  The 
Knudsen  number  is  defined  as  the  ratio  of  molecular  mean  free  path  to  a  characteristic 
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vehicle  dimension,  usually  nose  radius.    The  Knudsen  number  is  indexed  according  to 
its  reference  mean  free  path,  either  after  collision  or  in  the  oncoming  flow. 

The  boundary  layer  may  be  characterized  by  the  Reynolds  number 
(Re).  Flow  is  said  to  be  laminar  when  the  viscous  forces  are  sufficiently  large  to  damp 
out  oscillations  caused  by  the  dynamic  forces.  A  low  Reynolds  number  is  characteristic 
of  laminar  flow.  Turbulent  flow  is  said  to  occur  when  the  dynamic  forces  overcome  the 
viscous  forces,  there  is  random  mixing  of  particles  and  large  momentum  exchanges 
between  fluid  particles.  When  flow  over  a  solid  body  reaches  a  critical  Reynolds  number 
the  initial  laminar  flow  transitions  to  turbulent  flow.  Turbulent  flow  is  a  critical  factor 
in  reentry  since  there  is  much  more  energy  near  the  vehicle's  surface  than  in  laminar 
flow  conditions.  In  turn,  under  turbulent  flow  conditions,  more  heat  is  transferred  to  the 
surface  of  the  vehicle.  The  five  flow  regimes  are  shown  graphically  in  Figure  9  and 
described  as  follows:  [Ref.  20:pp.  203-206] 


1.  Free  Molecule  Regime-this  region  is  where  the  molecular  mean  free  path  (X) 
is  relatively  large  compared  to  a  characteristic  vehicle  dimension  such  as  nose 
radius.  When  molecules  collide  with  a  boundary  layer  they  attain  the  state  of 
that  boundary  after  a  single  collision.  This  flow  regime  is  characteristic  of  the 
uppermost  portion  of  the  atmosphere. 

2.  Near  Free  Molecular  Flow-frequently  referred  to  as  the  "slip  region"  because 
gas  molecules  will  acquire  the  momentum  of  the  moving  boundary  only  after 
several  collisions.  If,  on  the  average,  a  molecule  fails  to  acquire  the  momentum 
of  the  moving  boundary  after  a  single  collision,  then  it  is  said  to  lack 
"accommodation."  This  lack  of  accommodation  means  that  the  temperature  is 
a  nearly  discontinuous  function  of  distance  away  from  the  moving  surface. 
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Figure  9:    Flow  Regimes 
[Ref.  20] 
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3.  Transition  Regime—little  is  known  about  aerodynamic  quantities  such  as  lift, 
drag  or  heat  transfer.   This  region  is  not  treated  very  well  analytically. 

4.  Viscous  Merged  Layer  Regime— essentially  that  region  consisting  of  dynamically 
coupled  shock-wave,  boundary-layer  interactions.  The  presence  of  a  boundary 
layer  on  the  wall  alters  the  boundary  conditions  for  the  shock  wave, 
simultaneously,  the  large  pressure  gradients  due  to  the  shock  wave  strongly  alter 
the  boundary-layer  flow.  Neither  the  shock  wave  nor  the  boundary  layer  may 
be  treated  as  discontinuities.  This  is  the  region  from  approximately  1 10  km  to 
75  km;  this  is  where  the  "initial  pitch  over"  occurs  and  reentry  "starts." 

5.  Continuous  Regime-the  region  where  classical  fluid  mechanics  of  high  Reynolds 
number  applies,  here  the  shock  wave  and  boundary  layers  are  again  treated  as 
discontinuities.  Often  this  regime  is  subdivided  into  four  categories  (subsonic, 
transonic,  supersonic,  hypersonic)  with  lines  of  demarcation  established  by 
Mach  numbers. 


(4)  Vehicle  Profile  Area.  The  profile  area,  or  aerodynamically 
effective  cross-section  of  the  vehicle,  is  the  area  presented  to  the  oncoming  flow  of  the 
atmosphere  and  is  a  function  of  the  vehicle's  attitude  and  configuration.  A  spherical 
satellite  maintains  a  constant  profile  area.  More  complex  vehicles  that  have  various 
design  shapes  and  deployed  solar  panels  can  have  highly  variable  areas.  The  ability  to 
model  the  area  depends  on  the  known  dimensions  and  orientation  of  the  vehicle.  Since 
an  uncontrolled  satellite  reentry  is  characterized  by  a  critical  system  failure  or  the  end 
of  its  operational  lifetime,  communications  with  the  vehicle  may  have  been  lost.  Under 
this  condition,  information  on  the  satellite's  attitude  will  not  be  directly  available  [Ref. 
37:p.  4].  Additionally,  the  vehicle  will  experience  increasing  aerodynamic  and  thermal 
loads  as  the  altitude  decreases.  These  forces  will  act  to  change  the  vehicle  configuration 
by  deformation  and  removal  of  structures  (mass)  [Ref.  5:p.  iv]. 
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(5)  Vehicle  Mass.  Knowing  the  mass  of  a  satellite  accurately  assumes 
a  comprehensive  knowledge  of  the  vehicle.  This  information  may  be  available  for  some 
U.S.  satellites,  however,  for  foreign  systems  this  may  pose  a  problem.  As  mentioned 
in  the  previous  paragraph,  aerodynamic  and  thermal  loads  can  change  the  satellite  mass. 
The  ability  to  model  changes  in  mass  also  assumes  an  accurate  knowledge  of  the 
environment  and  how  it  affects  the  forces  and  subsequent  motion  of  the  vehicle  about  its 
center  of  mass.  In  order  to  develop  a  model  capable  of  predicting  this  motion,  the 
vehicle's  moments  of  inertia  must  be  known.  [Ref.  37:p.  4] 

(6)  Vehicle  Velocity.  The  relative  velocity  of  the  vehicle  with  respect 
to  the  Earth's  rotating  atmosphere  in  equation  (6)  is 

VA  =  \V-aXr\v  (11) 

The  maximum  deceleration  and  heating  rate  experienced  by  a  reentry  body  is  a  function 
of  velocity.  As  mentioned  before,  wind  can  change  the  vehicle's  relative  airspeed  which 
can  affect  the  drag  or  lift  experienced  by  the  vehicle.  Determination  of  the  vehicle's 
velocity  using  doppler  range  rate  observation  information  is  relatively  accurate. 
Uncertainties  in  the  radial  velocity  may  be  as  low  as  .166  m/s  [Ref.  38:p.  138]. 

d.    Modeling  Aerodynamic  Lift 

Lift  causes  the  vehicle  trajectory  to  follow  a  glide  slope,  skip  path,  or 
perturbed  ballistic  flight  path.  During  hypersonic  flight  conditions,  lift  is  generated  by 
pressure  forces  on  the  lower  surfaces  at  angles  of  attack  caused  by  motion  about  the 
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center  of  mass.  A  perturbation  from  the  nominal  or  zero  lift  flight  path  is  the  net  result. 
The  force  of  aerodynamic  lift  is  defined  as  [Ref.  39:p.  25] 

L  =  ipQAV2  (12) 

2      L 

where 

CL  =  coefficient  of  lift 
Reentry  trajectories  with  lift  reduce  the  thermal  and  structural  loads  on  the  vehicle. 
During  the  final  revolutions  of  orbit  decay  and  reentry,  lift  is  assumed 
or  modeled  to  be  very  small  or  zero  [Ref.  14:p.  2],  [Ref.  39:p.  33].  However,  an 
uncontrolled  reentering  satellite  may  generate  a  significant  lift  vector  depending  on  the 
attitude,  shape  and  motion  of  the  vehicle  about  its  center  of  mass.  When  this  occurs,  the 
lift  will  not  be  distributed  equally  about  the  nominal  flight  path.  This  may  result  in  a 
deviation  from  the  projected  impact  point.  [Ref.  15:p.  47] 

e.     Reentry  Equations  of  Motion 

The  reentry  equations  of  motion  can  be  derived  by  the  mathematical 

transformation  of  the  second-order,  nonlinear,  vector  differential  equation  (4).   Several 

analytical  and  semi-analytical  theories  describing  a  satellite's  shallow  reentry  equations 

of  motion  have  been  developed  over  the  last  thirty  five  years  [Refs.  14,40,41,42]. 

With  the  advent  of  manned  space  exploration  and  recoverable  probes,  it  became 
necessary  to  develop  accurate  theories  of  the  entry  phase,  during  which  the 
altitude,  velocity,  deceleration  and  the  heating  rate  vary  rapidly.  [Ref.  41:p.  2] 

In  the  derivations  of  these  theories,  several  strong  physical  assumptions  were  made: 
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1.  Spherical  Symmetry—the  Earth  and  its  atmosphere  are  spherically  symmetric. 

2.  Non-rotating  Atmosphere-the  rotation  of  the  Earth's  atmosphere  which  is 
approximately  equal  to  the  angular  velocity  of  the  planet  is  neglected. 

3.  Exponential  Atmosphere-the  atmospheric  density  decreases  exponentially  with 
altitude. 

4.  Gravitational  Field— the  gravitational  field  is  assumed  to  be  constant  during 
reentry  at  all  altitudes. 

5.  Coordinate  System— a  non-rotating  two-dimensional  inertial  coordinate  system 
with  the  origin  at  the  center  of  the  Earth. 

(1)  Fundamental  Equations  of  Entry  Dynamics.  The  exact  reentry 
equations  of  motion  derived  by  Loh,  reference  [40],  are  presented  to  show  the  basic 
relationship  of  the  vehicle's  force  vectors  along  the  radial  and  normal  direction  to  the 
flight  path  in  an  inertial  coordinate  system  as  shown  in  Figure  10  [Ref.  40:p.  19].  This 
presentation  serves  as  the  foundation  for  the  development  of  more  sophisticated  theories 
to  be  presented  in  Chapter  III.  Physical  assumptions  1,3,4,  and  5  were  used  in  the 
derivation  of  these  equations. 

The  magnitude  of  the  relative  aerodynamic  velocity, Vr ,  does  not 
equal  the  magnitude  of  the  inertial  velocity,  V,  due  to  the  reentry  body's  trajectory 
through  the  moving  atmosphere.  These  velocities  and  their  geometric  relationships  are 
shown  in  Figure  11   [Ref.  40:p.  18]  and  are  related  by  the  following  equation 
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Figure  10:  Inertial  Coordinate  System 
[Ref.  40] 
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Figure  11:  Relationship  Between  Relative  And  Inertial  Velocity 
[Ref.  40] 
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where 


Ve      =  1519  cos  5  (ft/ sec)  (velocity  of  the  Earth's  surface  at  a  specific  latitude) 

5        =  latitude  angle 

7        =  0  (flight  path  angle  to  the  local  horizon  from  Figure  10) 

The  force  components  from  Figure  10  are  defined  as 


F,  =  -Dcos  4>+Isin  ty-mgsm  iJt  =  m — (Kcos  4>) 


dt 


(14) 


F    =  Lcos  <J)  +  Dsin  ty-mgcos  i|r  =  m — (Ksin  4>)  (15) 

where 

4>       =  7  +  \p 

V       =  Vr  (relative  velocity) 
By  resolving  equations  (14)  and  (15)  in  the  velocity  direction,  replacing  L  and  D  with 
equations  (12)  and  (6)  respectively,  and  by  further  rearrangement  of  the  terms  and 
simplification,  according  to  Loh,  the  exact  equations  of  motion  are 
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where 

Ro      =  radius  of  the  Earth 

£        =  inverse  atmospheric  scale  height  (Mg/RT) 

M      =  mean  molecular  weight 

The  exact  equations  of  motion  of  entry  dynamics  cannot  be  solved 
analytically,  however,  solutions  can  be  found  using  numerical  methods.  First-order 
approximate  analytical  solutions  have  been  developed  by  restricting  the  equations  to 
limited  regions  of  application.  [Ref.  43:p.  25]  The  ability  to  accurately  model  and  solve 
these  equations  depends  on  the  knowledge  of  the  initial  conditions:  the  initial  position, 
and  velocity,  the  vehicle's  area  and  mass,  the  neutral  atmospheric  density,  and  the 
aerodynamic  coefficients. 

2.      Modeling  Breakup 

a.     Reentry  Body  Structural  Mechanics 

A  satellite  will  experience  structural  and  thermal  loads  during  reentry  into 
the  Earth's  atmosphere.  The  structural  response  of  the  body  during  reentry  may  depend 
on  the  following  factors:  [Ref.   44] 

1.  Static  load  effects 

2.  Dynamic  load  effects 

3.  Thermal  load  effects 
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The  coupled  effect  of  these  factors  may  determine  when  the  body  will  experience 
structural  failure. 

b.    Modeling  Reentry  Heating  Effects 

The  reentry  process  is  essentially  that  region  of  flight  where  the  vehicle's 
velocity  into  the  atmosphere  is  reduced  and  its  kinetic  energy  is  converted  into  thermal 
energy  in  the  surrounding  medium.  Since  breakup  is  determined,  in  large  part,  by  when 
the  outer  surface  reaches  its  melting  point,  reentry  heating  directly  affects  survivability. 
This  will  be  discussed  in  detail  in  Chapter  III.  The  conversion  fraction  of  kinetic  energy 
to  heat  energy  is  a  function  of  the  satellite's  shape,  velocity,  and  altitude.  At  very  high 
altitudes  (free-molecule  flow  region)  the  heat  energy  is  developed  almost  entirely  at  the 
surface  of  the  vehicle  and  up  to  one-half  of  the  lost  kinetic  energy  may  be  converted  into 
heat  in  the  vehicle  body.  [Ref.  43 :p.  191]  At  low  altitudes  (continuum  flow  region)  the 
heat  energy  will  appear  in  the  area  between  the  shock  wave  and  the  body.  This  heat  is 
transferred  from  the  hot  gas  to  the  vehicle  by  conduction  and  convection  through  the 
viscous  boundary  layer  which  is  adjacent  to  the  surface  of  the  vehicle.  Radiant  heating 
of  the  vehicle  from  the  hot  gas  also  occurs  and  this  contribution  to  the  surface  heating 
is  dependent  upon: 

1.  Bluntness  of  the  vehicle's  leading  edge 

2.  Excess  orbital  velocity  of  the  vehicle 
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Excess  orbital  velocity  is  defined  as  an  appreciably  higher  velocity  than  the  circular 
orbital  velocity  for  a  given  altitude.  [Ref.  43:p.  192] 

The  rate  at  which  the  vehicle  heats  is  not  exclusively  dependent  upon  this 
energy  conversion  fraction,  it  is  also  dependent  upon  the  rate  of  kinetic  energy  loss  by 
the  vehicle.  For  the  case  of  natural  orbital  decay,  reentry  is  at  small  flight-path  angles 
and  the  deceleration  is  very  slow  in  the  upper  atmosphere.  The  surface  heating  rate  is 
relatively  low  despite  the  high  conversion  fraction  in  this  instance;  therefore,  the 
dominating  factor  is  the  low  rate  of  kinetic  energy  loss  by  the  vehicle.  Conversely,  at 
steep  reentry  angles,  where  deceleration  occurs  rapidly,  the  surface  heating  rate  is  high 
although  the  energy  conversion  fraction  is  low. 

The  total  heat  input  into  the  reentry  vehicle  depends  upon  the  time  of  heating 
as  well  as  the  heating  rate.  If  the  energy  conversion  rate  were  constant,  the  total  heat 
input  would  simply  be  a  fixed  fraction  of  the  initial  energy  and  the  type  of  reentry  would 
not  be  of  significance. 

Three  aspects  of  the  aerodynamic  heating  process  are  significant,  namely: 
[Ref.  40:p.  181] 

1.  The  total  heat  input,  Q. 

2.  The  time  rate  and  maximum  time  rate  of  local  stagnation  region  heat  input  per 
unit  area  (dH/dt)  and  (dH/dt),^. 

3.  The  time  rate  and  maximum  time  rate  of  average  heat  input  per  unit  area 
(dH.v/dt)  and  (dJVdtU. 
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The  time  rate  of  average  heat  transfer  per  unit  area  is  given  by  equation  (18)  below.  The 
average  heat  transfer  is  simply  the  total  heat  input  divided  by  the  time  of  input.  Overall, 
reentry  vehicle  structural  integrity  is  a  function  of  the  average  heat  input.  [Ref.  40:p. 
181]  The  time  rate  of  local  stagnation  region  heat  input  per  unit  area  is  given  by 
equation  (19).  Local  structural  integrity  is  a  function  of  local  stagnation  region  heat 
input  or  the  generation  of  local  "hot  spots."  [Ref.  40:p.  182]  The  total  heat  input  is 
given  by  equation  (20). 

(1)  Reentry  Heat  Input.  From  the  three  previously  mentioned  areas  of 
concern,  the  following  generalized  aerodynamic  heating  equations  are  given  [Ref.  40:pp. 
183-184] 
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CF'  =  equivalent  skin-friction  coefficient 

V  =  velocity,  ft/sec 

p  =  atmospheric  density,  slugs/ft3 

Kn  =  Knudsen  number  =  X/Re 

Re  =  Reynolds  number 

X  =  molecular  mean  free  path 

H  =  coefficient  of  viscosity,  lb-sec/ft2 

Pr  =  Prandtl  number,  subscript  e  indicates  "entry"  value 

7  =  ratio  of  specific  heats 

a  =  nose  or  leading  edge  radius,  ft 

K'  =  constant  =  6.8  to  15  x  10^ 
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=  surface  area,  ft2 

=  flight  path  angle,  positive  for  descent 
=  coefficient  of  friction  (local  conditions) 
=  specific  heat  at  constant  pressure 


(22) 
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Cp,     =  specific  heat  (local  conditions) 
The  above  equations  make  the  following  simplifying  assumptions:  [Ref.  40:p.  181] 

1 .  Radiative  heat  transfer  from  the  surface  generally  does  not  appreciably  influence 
convective  heat  transfer  to  a  vehicle;  therefore,  it  is  disregarded. 

2.  Effects  of  gaseous  imperfections  may  be  neglected. 

3.  Shock- wave  boundary-layer  interactions  may  be  neglected. 

4.  Prandtl  number  is  constant. 

5.  Reynolds  analogy  is  applicable. 

If  the  heat  transferred  to  the  reentry  vehicle  is  expressed  as  a 
fraction  of  the  total  kinetic  energy,  then  [Ref.  20:p.  138] 

KE  =  ±mV2E  (23) 

2      E 

where 

m       =  mass  of  reentry  vehicle 
VE      =  reentry  velocity 

Q  =  -^  (24) 


where 

Q       =  fraction  of  heat  transferred  to  the  reentry  vehicle 
Q       =  total  heat  input 
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This  may  be  rewritten  as 


q  « 1  Eiz  [l-g-N 

*       2  C^S      L  J 

where 

CF  =  coefficient  of  friction 

CD  =  coefficient  of  drag 

S  =  reference  area  (usually  nose  tip) 

Sw  =  wetted  surface  area 

and  where  the  trajectory  parameter  is  defined  by  [Ref.  20:p.l38] 
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7E      =  reentry  flight  path  angle 
At  small  flight  path  angles,  such  as  reentry  from  orbit  decay,  aE  will  be  very  large, 
therefore,   equation  (25)  reduces  to 
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It  is  now  recognizable  that  the  fraction  of  the  initial  kinetic  energy,  transferred  to  the 
reentry  vehicle  by  convective  heating,  is  one  half  of  the  ratio  of  the  friction  drag  to  total 
drag  .  [Ref.  20:p.  138] 

In  the  case  of  ballistic  reentry  at  small  angles  of  reentry,  both  the 
maximum  heating  rate  and  total  heat  load  increase  as  the  effective  mass-area  ratio 
(m/CDA)  or  ballistic  coefficient  increases  [Ref.  40:p.  195]. 

Assuming  a  constant  ballistic  coefficient  throughout  the  reentry  process 
yields  a  different  relationship  [Ref.  40:p.  198] 
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Therefore,  reentry  at  small  angles  of  inclination  reduces  the  maximum  heating  rate  of  the 
vehicle,  however,  the  total  heat  load  is  increased.  [Ref.  40:p.  198] 
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(2)  Reentry  Heating  Rate.  If  the  heating  rate  per  unit  time  is  defined 
[Ref.  20:p.  135] 

d  =  *1  (32) 

dt 

then  the  total  heat-transfer  rate  may  be  written  [Ref.  20:p.  136] 

Q  =  /  qds  (33) 

jsw 

After  some  simplification  it  is  possible  to  write  [Ref.  20:p.  137] 

.  =  cF  P  sw  v3  (34) 

*  4 

this  assumes 

\sC/ts  =  C^w  (35) 

It  is  now  possible  to  define  dq/dt  as  the  average  rate  of  change  of  heat  transfer  per  unit 
area 

4  =  ^Z!  (36) 

q  4 

It  can  be  shown  that  dq/dt  is  a  maximum  when  pV3  is  a  maximum.  [Ref.  20:p.  137] 

(3)  Ablation.  As  previously  discussed,  the  heating  experienced  by 
reentry  at  small  flight  path  angles  is  different  from  that  of  reentry  at  large  flight  path 
angles.  The  flight  duration  is  much  longer  for  the  first  case  and  even  with  a  lower 
maximum  heating  rate,  the  total  heat  input  exceeds  that  of  the  larger  reentry  angle.  The 
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predominant  cooling  mechanism  of  uncontrolled  reentry  of  vehicles  not  designed  to 

survive  reentry  is  ablation.   Ablation  is  the  term  which  generally  applies  when: 

...there  is  a  removal  of  material  (and  an  associated  removal  of  heat)  caused  by 
aerodynamic  heating,  and  therefore  embraces,  melting,  sublimation,  melting  and 
subsequent  vaporization  of  the  liquid  film,  burning  and  depolymerization.  [Ref. 
15:p.  198] 

Ablative  materials  are  measured  as  "effective"  depending  on  their 
capacity  to  dispose  of  heat  (latent  heat)  by  convection  in  the  liquid  film,  and  by 
convection  in  gaseous  form  in  the  boundary  layer.  Previous  work  has  shown  that  when 
a  material  experiences  a  large  percentage  of  total  mass  loss  as  vaporization,  it  is  a  more 
effective  ablative  material.  Sublimation  is  the  process  whereby  all  the  mass  loss 
undergoes  vaporization,  there  is  no  liquid  film,  and  therefore  is  an  excellent  method  for 
removing  large  amounts  of  heat.  For  these  reasons,  materials  which  undergo  sublimation 
at  reasonably  high  temperatures  are  generally  more  efficient  at  removing  heat  and  thereby 
reducing  the  thermal  load  on  the  reentry  vehicle.  [Ref.  40:p.  204] 

Radiative  cooling  is  another  mechanism  by  which  heat  is  removed 
from  the  reentry  vehicle  during  reentry.  Radiative  cooling  and  ablation  combined  are  an 
effective  pair  in  balancing  the  heating  effects  of  reentry  for  a  lifting  body.  For  a  non- 
lifting  body  radiative  cooling  is  inefficient. 

C.      GENERAL  DESCRIPTION  OF  THE  REENTRY  PHASE 

A  decaying  satellite  possesses  a  large  amount  of  kinetic  energy  due  to  its  velocity 
and  potential  energy  because  of  its  altitude  above  the  Earth's  surface.  As  the  satellite 
encounters  the  atmosphere,  a  shock  wave  forms  ahead  of  the  vehicle  which  heats  up  the 
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atmosphere  surrounding  the  vehicle.  This  enveloping  layer  of  incandescent  atmosphere 
causes  the  vehicle's  temperature  or  thermal  load  to  continually  increase  as  it  penetrates 
into  an  increasingly  denser  atmosphere.  During  this  phase,  the  velocity  continues  to 
decrease  as  the  kinetic  energy  is  converted  into  heat  through  the  atmospheric  drag.  If 
all  the  satellite's  energy  were  converted  into  heat  and  contained  within  the  vehicle,  then 
there  would  be  more  than  enough  energy  to  vaporize  the  vehicle.  However,  this  is  not 
the  case.  A  large  part  of  the  total  energy  is  diverted  away  from  the  vehicle  by  two 
processes.  The  first  process  unloads  a  major  fraction  of  the  heat  into  the  atmosphere  by 
a  strong  shock  wave  mechanism.  The  second  process  involves  the  radiation  of  heat  away 
from  the  hot  surfaces.  [Ref.  43:pp.  1-2] 

The  vehicle  also  experiences  structural  loads,  which  are  a  combination  of 
aerodynamic  and  inertial  loads  [Ref.  45 :p.  5].  Atmospheric  drag  forces  usually  cause 
a  reduction  in  the  vehicle's  kinetic  energy.  Centrifugal  and  lift  forces  cause  accelerations 
normal  to  the  direction  of  the  motion. 

Aerodynamic  lift  and  drag  forces  vary  directly  with  the  square  of  the  vehicle's 
velocity,  V2,  and  with  the  atmospheric  density,  p.  Deceleration  of  the  vehicle  is  a 
product  of  two  quantities,  density  and  velocity.  As  the  satellite  penetrates  further  into 
the  atmosphere,  the  density  increases  rapidly  resulting  in  a  corresponding  decrease  in 
velocity  due  to  drag.  Initially,  deceleration  increases  as  shown  in  Figure  12  [Ref.  43:p. 
7].  However,  at  some  altitude  the  velocity  begins  to  decrease  at  a  faster  rate  than  the 
increasing  density,  which  results  in  a  maximum  deceleration.  Additionally,  the 
maximum  heating  rate,  -  pV3,  occurs  at  an  altitude  somewhat  higher  than  the  maximum 
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Figure  12:  Changes  During  Atmospheric  Reentry 
[Ref.  43] 
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deceleration  [Ref.  20:p.  98],  [Ref.  43:pp.  6-7].  A  maximum  deceleration  of  up  to  8 
g's,  as  shown  in  Figure  13,  [Ref.  14:p.  30],  occurs  during  this  phase  [Ref.  39:p.  36]. 
When  the  satellite's  skin  temperature  and  structural  load  become  sufficiently  high, 
the  vehicle  will  start  burning  and  breaking  up.  Solar  panels  and  other  projections  such 
as  antennas  will  separate  at  the  earliest  stage,  while  heavier  pieces  will  breakup  at  lower 
altitudes.  Based  on  predicted  and  observed  data,  satellite  breakup  commences  at  an 
altitude  between  approximately  75-120  km  [Ref.  10:p.  5],  [Ref.  46:p.  4],  [Ref.  47:p. 
39].  The  resulting  debris  from  the  breakup  will  impact  the  Earth's  surface  provided  it 
survives  the  reentry  heating  process. 

D.      CURRENT  REENTRY  THEORIES 

Predicting  reentry  time  and  impact  location  relies  upon  observational  data  of  the 
reentry  vehicle.  Based  upon  observed  (measured)  position  and  velocity,  over  the  orbital 
path  and  ideally  equally  distributed  over  that  path,  an  algorithm  is  used  to  calculate  an 
elliptical  orbit  which  best  fits  the  observational  data.  If  the  algorithm  attempts  to  model 
the  physical  reentry  process,  it  will  be  referred  to  as  a  physical  model  by  the  authors. 
Another  type  of  model  to  be  discussed  is  the  King-Hele  or  mean  motion  type  of  model 
which  neglects  certain  physical  aspects  of  the  reentry  process  and  focuses  more  on  the 
observational  data.  The  algorithm  which  "fits"  the  observational  data  to  an  elliptical 
orbit  may  also  be  used  to  propagate  or  predict  future  orbital  locations  as  a  function  of 
time.   These  types  of  algorithms  are  commonly  referred  to  as  propagators. 
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1.       Physical  Modeling  Theory 

The  current  propagator  in  the  United  States,  for  reentry  prediction,  is  the 
special  perturbations  (SP)  model.  This  model  is  maintained  by  the  Air  Force  Space 
Command  and  is  the  standard  for  reentry  prediction  in  the  United  States. 

Orbital  periods  of  less  than  87.5  minutes  are  defined  by  the  Space 
Surveillance  Center  (SSC),  Cheyenne  Mountain  Air  Force  Base,  Colorado,  as  decaying 
orbits  [Ref.  48:p.  3].  The  SP  model  uses  numerical  methods  to  incorporate  zonal, 
sectoral  and  tesseral  orbital  perturbations  in  the  calculation  of  decaying  orbit  reentry 
predictions.  Gravitational  perturbations  are  modeled  by  mapping  the  Earth  into  small 
grids,  which  allows  for  enhanced  resolution  of  the  geopotential.  Third  body  gravitational 
effects,  sun,  moon,  and  planets,  may  also  be  modeled.  Third  body  gravitational  effects 
are  used  primarily  in  the  propagation  of  highly  eccentric  orbits  with  large  apogee 
altitudes.  The  atmosphere  is  modeled  using  the  Jacchia-Nicolet  model  which  considers 
the  following:  [Ref.  48:p.  9] 

1.  Diurnal  bulge 

2.  Solar  activity 

3.  Geomagnetic  activity 

4.  Semiannual  variation 

A  technique  known  as  differential  corrections  is  used  to  "fit"  the  observations 
of  the  reentry  vehicle  into  the  best  orbit.  The  differential  corrections  are  a  mathematical 
means  of  determining  a  single  orbit  path  (ellipse)  consistent  with  the  observed  data  (the 
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only  "known"  information  in  the  absence  of  active  telemetry  data).  The  need  for  the  use 
of  differential  corrections  arises  from  the  fact  that  the  model  used  to  propagate  the  orbit 
into  the  future  has  inherent  deficiencies.  These  deficiencies  are  then  reflected  as  errors 
in  the  prediction  of  where  the  vehicle  is  supposed  to  be  in  its  orbit  as  compared  to  where 
it  is  observed  to  be.  Through  a  series  of  fitting  observations  to  an  orbit  and  updating  the 
predicted  orbital  path,  the  reentry  of  the  vehicle  is  calculated  as  a  "time"  when  the 
altitude  will  reach  a  specified  minimum  value.  This  lower  limit  value  of  altitude  is  a 
function  of  the  atmospheric  density  model.  The  impact  dispersion  area  or  footprint  is 
calculated  as  ±  15  minutes  of  that  reentry  time.  The  sub-satellite  ground  trace  is  used 
to  describe  the  impact  area  on  the  Earth's  surface.  Table  2  shows  the  format  of 
observational  data  as  used  in  the  differential  corrections  process  [Ref.  49].  Figure  14 
shows  the  differential  corrections  display  and  Figure  15  shows  the  Tracking  and  Impact 
Prediction  (TIP)  display  used  at  Cheyenne  Mountain  AFB.  [Ref.  49] 

It  now  becomes  obvious  that  there  are  two  extremely  significant  issues  at 
hand: 


1.  The   ±15  minute  window  equates  to  approximately   1/3  of  one  complete 
revolution  of  the  reentry  vehicle's  orbit  in  the  decay  phase. 

2.  The  lack  of  observational  data  in  the  decay  phase  may  significantly  bias  the 
predicted  reentry  time. 
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Table  II: 

[Ref.  49] 
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Figure  14:  Differential  Corrections  Display 
[Ref.  49] 
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Figure  15:  Tracking  And  Impact  Prediction  (TIP)  Display 
[Ref.  49] 
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The  accuracy  of  the  reentry  prediction  is  directly  limited  by  the  ability  to 
observe  the  reentry  process  as  well  as  indirectly  by  the  inherent  deficiencies  in  the 
models  used  to  represent  the  physical  reentry  process. 

2.       Mean  Motion  Theory 

In  the  early  development  of  his  mean  motion  theory  of  satellite  lifetime 
prediction,  King-Hele  defined  orbital  lifetime  as  the  time  remaining  until  the  eccentricity 
of  the  orbit  reached  zero  [Refs.  19,  50-55].  This  was  a  prediction  of  reentry  based  on 
the  circularization,  or  contraction  of  the  orbit,  due  to  atmospheric  drag.  A  principal 
assumption  in  this  theory  was  that  perigee  height  remained  constant.  When  oscillations 
in  perigee  height  were  considered  (which  occurs  when  an  oblate  atmosphere  is  modeled), 
the  theory  had  to  be  revised  in  order  to  take  into  account  a  zero  eccentricity  prior  to 
reentry  [Ref.  54].  The  definition  of  orbital  lifetime  was  then  redefined  in  terms  of  mean 
motion,  n,  where  a  value  of  n  (chosen  by  the  user)  represents  an  "end  value"  of  a 
circular  height.  For  example,  a  value  of  n  equal  to  16.5  rev/day  would  correspond  to 
a  circular  orbit  altitude  of  150  km. 

In  the  prediction  of  orbital  lifetime,  the  observed  rate  of  change  of  orbital 
period,  dT/dt,  was  the  entering  argument  (where  T  is  the  orbital  period)  [Refs.  19,  50- 
53].  This  was  equivalent  to  using  the  rate  of  change  of  mean  motion,  dn/dt,  since  the 
two  are  related  by  [Ref.  54 :p.  5] 


n  =  I  (37) 

T 
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where 

n        =  rev/day 

T        =  days 
In  its  simplest  form,  orbital  lifetime  is  now  expressed  in  terms  of  the  rate  of  change  of 
mean  motion  as 

L  =  5  (38) 

n 

where 

L        =  lifetime  in  days 

dn/dt  =  rev/day2 

Q       =  function  of  e,  n  and  H  (density  scale  height) 
and  expanded  by  [Ref.  54 :p.  8] 


e=^[l-0.17]  (39) 

47,(2)  L  J 


where 

IojIj    =  Bessel  functions  of  the  first  kind,  of  degree  0  and  1,  with  argument  z 
J        =  0.3eo  -  0.025 
z        =  ae/H 

For  the  circular  orbit,  the  lifetime  may  be  written  [Ref.  54:p.  9] 
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lafi 
where 

H_H     =  scale  height  H  km  below  initial  altitude 
h        =  decrease  in  altitude  until  reentry  (120  km) 

In  this  form,  density  scale  height  is  the  only  parameter  which  is  not  directly 
derived  from  orbital  data  (observations).  It  can  be  shown  that  H  is  less  sensitive  to  solar 
activity  than  a  neutral  density  atmosphere  model,  therefore,  uncertainties  in  solar  activity 
will  have  less  influence  on  the  orbital  lifetime  prediction  [Ref.  54:pp.  9-10].  Perhaps 
more  importantly,  the  mean  motion  technique  does  not  require  any  knowledge  of: 

1.  Satellite  attitude 

2.  Satellite  size 

3.  Satellite  shape 

4.  Satellite  mass 

5.  Aerodynamic  characteristics 


These  factors  are  all  accounted  for  in  the  observation  of  mean  motion  and  rate 
of  change  of  mean  motion.  In  summary,  the  mean  motion  theory  is  dependent  upon  five 
parameters:  [Ref.  54:p.  12] 
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1.  mean  motion  (n) 

2.  eccentricity  (e) 

3.  daily  solar  activity  index  (F,07) 

4.  average  solar  activity  index  (F107) 

5.  geomagnetic  index  (A,,) 


The  last  three  terms  above  are  incorporated  into  the  calculation  of  atmospheric  scale 
height. 
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m.   FORMULATIONS  AND  SOLUTIONS  OF  REENTRY 

This  chapter  describes  the  "state-of-the-art"  formulations  and  solutions  (analytical, 
semi-analytical,  and  numerical)  of  decay-induced  reentry  as  determined  from  the 
literature  survey.  Fundamental  physical  processes  that  were  introduced  in  Chapter  II, 
including  reentry  equations  of  motion,  rarefied  gas  dynamics,  reentry  heating  and  reentry 
body  breakup  are  presented. 

A.      ANALYTICAL  REENTRY  EQUATIONS  OF  MOTION 

1.       Chapman's  Approximate  Analytical  Entry  Equations  of  Motion 

Chapman  derived  a  nonlinear  second-order  differential  equation  by  reducing 
equations  (16)  and  (17)  and  by  introducing  a  set  of  completely  nondimensionalized 
variables  [Ref.  14:pp.  3-14],  In  the  development  of  this  equation,  physical  assumptions 
1  through  3  from  Chapter  II  were  used.  Additionally,  two  mathematical  assumptions 
were  used: 


1.  The  fractional  change  from  the  center  of  the  planet,  dr/r,   in  a  given  increment 
of  time  is  much  smaller  as  compared  to  the  fractional  change  in   velocity,  d(V 
cos  7)/V  COS7,  given  by  the  mathematical  expression    |  d(V  cos  7)/V  cos  7  | 
>    I  dr/r  I  . 

2.  The  flight  path  angle,  7,  related  to  the  local  horizon  direction  for  lifting  vehicles 
is  sufficiently  small  such  that  the  lift  component  in  the  horizontal  direction  is 
small  as  compared  to  the  drag  component,  given  by  the  mathematical 
expression,  1  >    |  (L/D)  tan  7  |  . 
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It  is  erroneous  to  think  that  these  assumptions  will  restrict  Chapman's  analysis  to 
entry  trajectories  with  small  flight  path  angles  and  small  lift-to-drag  ratios  as  many 
authors  have  believed.  On  the  contrary,  these  assumptions  applied  simultaneously, 
constitute  a  well  balanced  set  of  hypotheses  and  make  Chapman's  theory  applicable 
to  a  large  family  of  entry  trajectories.  [Ref.  42 :p.  179] 

As  previously  mentioned,   Chapman  introduced  a  set  of  dimensionless 

(independent  and  dependent)  variables.   The  independent  variable  is 


(41) 


V  cos 
u  = 

7        u 

yfgr~        u< 

and  the  dependent  variable  is 

z-  fiC^ 

2m 

r   - 
—  u 

(42) 


Chapman's  final  derivation  is  given  by 


d2Z 
dU2 


dZ 
du 


Z] 

u 


<'  ;  "2>  cos'  y 
uZ 


J@r    —  cos3  7 
D 


(43) 


1. 


where  /3  =  Mg/RT  and  the  numbers  associated  with  the  brackets  in  equation  (43) 
represent  the  following  physical  quantities: 

1.  Vertical  acceleration. 

2.  Vertical  component  of  the  drag  force. 

3.  Force  of  gravity  minus  the  centrifugal  force. 

4.  Lift  force. 
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For  the  specific  case  of  the  shallow  reentry  of  a  satellite  from  a  decaying 
orbit  where  the  flight  path  angle  is  very  small,  the  Z  function  equation  (43),  can  be 
written  in  the  following  form 


-d 
du 


1  dZ      Z" 


-  ILfii  +  M-r  A  =  0  (44) 

uZ  D 


Kdu       Uj 

with  the  initial  conditions  of 

ut  =  1    ,    % ,  -  0  (45) 

and  where  the  corresponding  boundary  conditions  for  decaying  orbits  are 

Z(l)  =  0    ,    Z'(l)  =  0  <46> 

Integrating  equation  (44)  by  numerical  techniques  for  each  designated  initial 
velocity,  initial  flight  path  angle  and  L/D  ratio  generates  a  solution.  The  results  of  the 
solution  are  applicable  to  any  vehicle  of  arbitrary  dimension,  size,  or  mass.  Figure  16 
shows  the  solution  graph  for  a  non-lifting  vehicle  (L/D=0)  [Ref.  14].  Figure  17  shows 
the  solution  graph  of  various  L/D  values  using  the  parameter  [Ref.  14] 


D 


(47) 


The  ratio  u/uc  for  the  ordinate  from  Figures  16  and  17  is  equal  to  equation  (41). 

Additionally,  the  following  engineering  quantities  of  interest  during  the 
reentry  process  can  be  calculated  from  the  Z  function: 
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Dimensionless  velocity,  u-—- 


Figure  16:  Z  Function  Solution  Graph  For  L/D=0 
[Ref.  14] 


{^r)uZ 


Dimensionless  velocity,  u-^j- 


Figure  17:  Z  Function  Solution  Graph  For  Various  L/D's 
[Ref.  14] 
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1 .  Deceleration 

2.  Descent  angle 

3.  Range 

4.  Time 

5.  Density- velocity  relationship 

6.  Dynamic  pressure 

7.  Reynolds  number 

8.  Heating  rate 

9.  Total  heat  absorbed 

2.       Loh's  Second  Order  Unified  Solution  of  Entry  Dynamics 

Loh  derived  a  general  second-order  solution  of  reentry  mechanics  that  covers 
the  entire  range  of  initial  flight  path  angles  and  L/D  ratios  [Ref.  40:pp.  25-36].  Figure 
18  shows  the  entire  range  of  the  second-order  solution  as  compared  with  several  other 
analytical  approximate  solutions  available  in  1963,  including  Chapman's  approximate 
analytical  theory  [Ref.  40: p.  26]. 

The  second-order  unified  solution  can  be  derived  from  equations  (16)  and 
(17)  by  using  approximations,  a  binomial  series  expansion,  and  integration  techniques 
for  a  small  flight  path  angle,7,  and  is  given  by  the  following  equation 


71 


E scope 


w 


+fc 


Figure  18:  Loh's  Second-Order  Range  Of  Applicability 
[Ref.  40]  y 
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For  the  specific  case  when  L/D=0,  equation  (48)  can  be  simplified  and 


rewritten  as 
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Figure  19  shows  a  comparison  between  Chapman's  first-order  solution  and 

Loh's  second-order  solution  for  a  small  L/D=0.1  [Ref.  40:p.  63].  Both  solutions  offer 

about  the  same  degree  of  accuracy. 

However,  the  first  order  solution  is  limited  by  the  conditions  that  (L/D)  tan  7  must 
be  smaller  than  1  and  the  initial  angle  of  inclination  yf  must  be  very  small;  the 
second  order  solution  does  not  suffer  these  particular  limitations.  When  L/D=0 
and  7f  a  0,  where  the  first  order  solution  was  not  available  previously,  the  second 
order  solution  offers  for  the  first  time  a  satisfactory  solution  in  this  region. [Ref. 
40:pp.  50-51] 

3.       Yaroshevskii's  Entry  Theory 

Yaroshevskii  developed  a  semi-analytical  reentry  trajectory  theory  which  was 

originally  published  in  a  Soviet  journal,  Kosmicheskie  Issledovaniya,  in  1964  [Ref.  42:pp. 

158-176].  [An  English  translation  of  this  article  was  not  available  to  the  authors.] 

Using  some  simplifying  assumptions,  he  derived  a  nonlinear  second-order 
differential  equation  which  can  be  integrated  analytically  by  using  series 
expansions.  To  some  extent,  Yaroshevskii' s  theory  is  a  special  case  of  a  more 
sophisticated  theory  developed  by  Chapman.  [Ref.  42:p.  158] 

In  the  development  of  the  theory,  physical  assumptions  1  through  5  apply  to 

the  basic  reentry  equations  (16)  and  (17),  from  which  the  differential  equation  was 

derived.    Also,  for  a  constant  angle  of  attack,  a,  the  drag  coefficient,  CD,  and  the  lift 

coefficient,  CL,  were  assumed  to  be  a  function  of  the  Mach  number. 
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Figure  19:  Loh's  Second-Order  vs  Chapman's  First-Order 
[Ref.  40] 
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Yaroshevskii  used  an  independent  variable,  x,  and  a  dependent  variable,  y, 
defined  as 


Ji     CD(V)V 
and 


CD(l)A 

y  s 


rQ  (54) 

1P 


2m 

where 

r0       =   radius  of  the  planet 


*  -  —  (55) 


2ro 

By  defining  the  differential  relations  dx  and  dy,  using  equations  (53)  and  (54)  along  with 
the  equations  (16)  and  (17)  and  by  eliminating  the  flight  path  angle,  7,  the  second-order, 
nonlinear  differential  equation  is  given  by 


1     -1 


cfy  _      r^-CL[V(x)]  V2(x)  (56) 

dx2         V     °      CD{\)  v 

Solutions  for  equation  (56)  can  be  obtained  by  numerical  integration.  When 

CD  and  CL  are  independent  of  the  Mach  number,  solutions  can  be  obtained  by  integrating 

a  selected  series  expansion  of  equation  (54),  depending  on  the  type  of  reentry  trajectory. 
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For  the  specific  case  of  the  shallow  reentry   from  orbit  decay,   where 
CL/CD=0,  equation  (56)  can  be  derived  as 

£>  -  £H  (57) 

dx2         y 
with  the  initial  conditions  of 

xt  -  0    ,    y(0)  =0,^=0  (58) 

dx 

By  taking  into  account  a  singularity  at  y=0,  the  series  solution  to  equation  (57)  is 

y  =     -  x1(a0  +  alx+a2x2+a3x:i  +...)  C59) 

where 

ao       =  1 

a,       =  1/6 

a2       =  1/24 

a3       =  47/4752 
The  coefficients  \  can  be  calculated  by  the  recurrence  formula 


1  t_1 
-l£(2m  +  l)(2>72+3)amfl 


at  = 


(*  +  !)!     3^  (60) 

1+^+lK2^3) 
3 


As  in  the  case  of  Chapman's  theory,  several  engineering  quantities  of  interest 
during  the  reentry  process  can  be  calculated  from  equation  (56): 
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1.  Time  of  flight 

2.  Range 

3.  Deceleration 

4.  Heating  rate 

5.  Heat  absorbed 

4.        Universal  Equations  for  Orbit  Decay  and  Reentry 

Longuski  and  Vinh  derived  a  set  of  universal  entry  equations  of  motion  for 
all  regimes  of  atmospheric  flight:  from  the  free-molecule  flow  regime  to  the  near  free- 
molecule  flow  regime  where  orbital  motion  is  perturbed  by  air  drag,  through  the 
transition  regime  to  the  continuum  flow  regime  where  the  dynamic  phase  of  reentry 
occurs,  and  to  the  point  of  impact  on  the  planet's  surface  [Ref.  41]. 

Rigorous  mathematical  techniques,  such  as  Poincare"s  method  of  small 
parameters,  and  Lagrange's  expansion  are  applied  to  obtain  a  highly  accurate, 
purely  analytical  theory  for  orbit  contraction  and  ballistic  entry  into  planetary 
atmospheres.  [Ref.  41:p.v] 

Figures  20  and  21  [Ref.  41:pp.  16-17]  describe  the  inertial  coordinate  system, 

nomenclature  and  the  aerodynamic  forces  of  the  equations  of  motion  for  a  vehicle  with 

a  CL/CD  ratio,  defined  by  the  following  equations  [Ref.  41:p.l0] 

^  =  Vsin7  (61) 

dt 


dd       Vcos  7  cos  \p 
dt  rcos  <f> 


(62) 
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Figure  20:  Coordinate  System  And  Nomenclature 
[Ref.  41] 
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Figure  21:  Aerodynamic  Force  Diagram 
[Ref.  41] 
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d<f>        Vcos  7  sin  \p 

dt  r 


(63) 


dV  P*CDV2 


dt 


2m 


gsm  y 


(64) 


dy        pAC^cos  a  V2 

V—    =   -(g- )CX>S  Y 

dt  2m  r 


(65) 


M.       pACLV2cos  a     v2  .  t      j 

V—L  =  1 —  COS7  cos\p  tan0 

dt  2/wcos  7         r 


(66) 


The  exact  universal  equations  of  motion  for  entry  trajectories  for  a  vehicle 
inside  a  rotating  atmosphere  can  be  derived  from  a  transformation  of  equations  (61) 
through  (66)  using  the  modified  Chapman  variables 


u  = 


^cos2 
gr 


(67) 


z  =  pAC° 

2m 


r 
1 


(68) 


and  nondimensional  independent  variable 


■j: 


V 
r 


cos  7  dt 


(69) 


The  exact  universal  equations  for  entry  trajectories  are  [Ref.  41:pp.  11-12] 
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(71) 


(72) 


da       ,  _  \firZ  sin  a 


<fr 


tan/  cos/7 


sin  a 


(73) 


tffl       ^gTz  sina 
ds         sin/  cos2 7 


sin  a 


(74) 


di_ 
ds 


\/pr    Z  cosa 


cos2  7 


cL 


sin  a 


(75) 


where 


« 


a        = 


l        = 


longitude  of  the  ascending  node  of  the  osculating  plane 
angle  between  the  ascending  node  and  the  position  vector 
inclination  of  the  orbit 


a       =  bank  angle 
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According  to  the  authors,  the  universal  equations  have  three  advantages:  [Ref. 
41:p.l28] 

1.  They  are  exact. 

2.  They  are  free  of  any  restrictive  assumptions. 

3.  They  contain  the  modified  Chapman  Z  variable,   which  permits  a  single 
trajectory  solution  for  a  specified  initial  velocity  and  flight  path  angle  that 
applies  to  any  vehicle  of  arbitrary  area,  mass,  or  CD 

As  previously  mentioned,  equations  (70)  through  (75)  can  be  used  during  all 
phases  of  aerodynamic  flight  and  are  most  useful  for  analyzing  the  last  few  revolutions 
and  the  reentry  phase.  The  accuracy  of  the  equations  depends  on  the  readjustment  of 
the  value  of  the  inverse  atmospheric  scale  height,  /3,  for  each  layer  of  the  atmosphere. 
[Ref.  41:p.l5] 

For  the  specific  case  of  reentry  from  a  circular  orbit,  where  CL=0,  Longuski 

and  Vinh  derived  a  separate  analytical  theory  from  the  exact  universal  equations  for  entry 

trajectories,  due  to  the  fact  that: 

...it  does  not  seem  possible  to  have  a  single  analytic  solution  which  is  uniformly 
valid  for  all  values  of  initial  flight  path  angles  because  of  the  nature  of  the 
problem.  In  the  case  of  atmospheric  entry  from  circular  orbit,  the  magnitude  of 
the  flight  path  angle,  initially  zero,  rapidly  increases,  approaching  90°  as  the 
velocity  becomes  small.  On  the  other  hand,  for  steep  angle  entry,  the  flight  path 
angle  changes  very  little—of  the  order  of  tenths  of  a  degree—as  the  nondimensional 
velocity  decreases  from  unity  to  one  tenth  of  the  original  value  (between  Mach  2 
and  3).  [Ref.  41:p.  85] 

Under  the  condition  of  reentry  from  a  circular  orbit,  equation  (73)  is  equal 
to  one  and  equations  (74)  and  (75)  are  equal  to  zero.  By  using  this  fact  and  the  variable 
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V   = 


(76) 


gr       cos2  7 


equations  (70)  through  (73)  can  be  written  as 


—  =  -0rZ  tan7 
da 


(77) 


dv_  m  _2y//3r~Zv 
da  cos 7 


±L  =  l-i 

da  v 


sin  7 


i-2 


(78) 


(79) 


By  dividing  equations  (77)  and  (79)  by  equation  (78)  to  form  a  new  set  of  equations,  and 
defining  the  following  change  of  variables  for  substitution  into  the  new  set 


Y  =  2Z 


(80) 


$  =  -ypr    sin7 
X  =  -In  v 


(81) 
(82) 


€    = 


0r 


(83) 


and  then  by  expanding  this  new  set  of  equations  in  one  term,  in  e,  and  using  Poincare"s 
method  of  small  parameters  for  integration  of  a  system  with  the  assumed  solution  form 
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<|)    =   <J>0      +      €<!>! 


(84) 


Y=Y0+eYl 


(85) 


two  systems  of  two  first-order  differential  equations  are  formed  [Ref.  41  :p.  89] 


dYt 
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dX 
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d*o  m  ex-l 

dX    "      Yn 


(87) 


dX 


=  $.+  —  (2ex-l) 


(88) 


d*x 
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^x"1 


-W-l)-*?-£ 


(89) 


Equations  (86)  and  (87)  can  written  in  the  following  second-order  differential  equation 
form 


d%  _  e*-l 

dX2    "      Yn 


(90) 


and  with  the  initial  conditions  for  the  case  of  the  shallow  satellite  reentry 


y0(0)  =  o  ,  #0(0)  =  o 


(91) 


a  series  solution  can  be  obtained  in  the  following  form 
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(93) 


Yaroshevskii's  approach  was  used  to  find  the  first  term  of  the  series  in  equations  (92) 
and  (93).  The  same  approach  can  be  applied  to  equations  (88)  and  (89)  to  find  a 
solution. 

Figures  22,  23,  and  24  represent  various  solutions  graphs  for  equations  (92) 
and  (93).  The  dashed  line  indicates  the  exact  numerical  solution  while  the  solid  line 
represents  the  analytical  solution.  Figure  22  shows  the  variations  of  aerodynamic 
deceleration,  G's  (g's  =  number  x  gravitational  acceleration  at  a  radial  distance  r),  as 
a  function  of  the  dimensionless  velocity,  v,  at  several  initial  flight  path  angles,^,  for 
equation  (92)  [Ref.  41  :p.  120].  Figure  23  shows  the  variations  of  In  (  Z  /  Zq)  « ( r0  -  r)/H 
=  drop  in  altitude,  in  units  of  scale  height,  as  a  function  of  (v)  at  several  (yt  ),  for 
equation  (92)  [Ref.  41:p.  122].  Figure  24  shows  the  variation  of  the  negative  flight  path 
angle,  -y,  as  a  function  of  (v)  at  several  (y{ ),  for  equation  (93)  [Ref.  41  :p.  121]. 

5.       Attitude  Dynamics  of  Uncontrolled  Motion  During  Reentry 

In  the  previous  section  several  analytical  theories  were  presented  describing 
the  reentry  equations  of  motion  and  their  solutions.  Strong  physical  assumptions  were 
made  in  the  derivations  of  these  theories  in  order  to  describe  the  trajectory  of  the  body's 
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Figure  22:  Variations  Of  -7  vs  The  Nondimensional  Velocity  (v) 
[Ref.  41] 
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Figure  23:  Variations  Of  ln(Z/Zo)  vs  The  Nondimensional  Velocity  (v) 
[Ref.  41] 
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Figure  24:  Variations  Of  G  vs  The  Nondimensional  Velocity  (v) 
[Ref.  41] 
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center  of  mass  or  point  mass.  Specifically,  Chapman's  and  Longuski  and  Vinh's  theories 

used  the  Z  variable  which  permits  a  single  trajectory  solution  for  a  specified  initial 

velocity  and  flight  path  angle  that  applies  to  any  vehicle  of  arbitrary  area,  mass,  or  CD. 

The  effect  of  the  uncontrolled  motion  of  a  body  about  its  center  of  mass  on  the  reentry 

trajectory  was  not  investigated  in  these  theories.     This  section  will  examine  three 

analytical  investigations  presented  by  two  Russian  authors  in  this  specific  area. 

Duzmak,  in  1970,  presented  the  first  systematized  explanation  on  the  problem 

of  the  attitude  dynamics  of  uncontrolled  reentry  body  motion.  This  problem  was  the  least 

developed  in  comparison  to  center  of  mass  or  point  mass  trajectories  and  aerodynamic 

heating  [Ref.  56:p.2] 

Unifying  papers  on  the  dynamics  of  the  motion  relative  to  the  center  of  mass  have 
not  appeared  up  to  the  present  time.  [Ref.  56:p.  5] 

The  primary  focus  was  the  investigation  and  derivation  of  the  relationship 

between  a  reentry  body's  parameters  outside  the  atmosphere  with  the  body's  parameters 

in  the  dense  layer  of  the  atmosphere.    Changes  in  the  state  of  a  reentry  body's  motion 

relative  to  its  center  of  mass  during  the  motion  along  its  trajectory  were  also  investigated. 

An  asymptotic  approach  was  used  on  the  equations  of  motion  to  solve  the  problem, 

specifically  for  the  cases  where  the  characteristic  time  of  the  entry  body  motion  relative 

to  its  center  of  mass  is  much  less  than  the  characteristic  time  of  motion  of  the  center  of 

mass.     Additionally,  a  refined  asymptotic  method  that  has  a  significantly  wider  range 

of  applicability  was  developed  for  those  cases  in  which  the  above  condition  was  not 

fulfilled.      This  method  is  based  on  the  coupling  of  the  numerical  and  asymptotic 
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solutions  of  the  equations  of  motion.  Two-dimensional  motion  without  restrictions  on 
the  shape  of  the  body  and  three-dimensional  motion  of  a  body  that  possessed 
aerodynamic  and  dynamic  axial  symmetry  were  assumed  in  this  investigation.  [Ref. 
56:pp.  1-3] 

During  the  orbital  or  exoatmospheric  phase,  the  external  moments  that 
produce  the  reentry  body's  perturbed  motion  about  its  center  of  mass  are  determined  by 
the  laws  of  motion  of  a  rigid  body  as  described  by  Euler.  The  magnitude  and  direction 
of  the  initial  angular  moment  vector,  H,  determines  the  state  or  nature  of  this  motion. 
For  instance,  if  initial  angular  moment  vector  during  the  orbital  phase  produces  two- 
dimensional  motion,  then  the  body  rotates  at  a  constant  angular  velocity,  w,  about  its  Z 
or  transverse  axis.  [Ref.  56:pp.  5-6]  Figure  25  is  an  attitude  dynamics  coordinate  system 
that  shows  the  direction  of  the  Z-axis  relative  to  the  body's  orbit  [Ref.  57:p.  113]. 

Regardless  of  the  nature  of  the  motion,  the  entry  body's  angle  of  attack,  a,  upon 
approaching  the  atmosphere  can  have  absolutely  any  value  at  all  ...  Thus  any 
complete  solution  of  the  problem  of  atmospheric  entry  can  be  obtained  upon  the 
necessary  condition  of  the  absence  of  any  restrictions  on  the  size  of  the  angle  of 
attack.  One  can  say  that  the  presence  of  large  angle  of  attack  is  one  of  the  main 
distinctive  features  of  the  problem  of  atmospheric  entry.  [Ref.  56:p.  6] 

As  a  body  with  perturbed  motion  relative  to  its  center  of  mass  enters  the 

reentry  or  atmospheric  phase,  it  begins  to  experience  a  stabilizing  effect  that  is 

proportional  to  the  increase  in  atmospheric  density.  This  stabilizing  effect  is  a  property 

of  dynamical  systems  with  variable  parameters  where  if  the  stiffness  of  the  system 

increases  during  the  perturbed  motion,  then  the  vibrations  dampen  out.  This  dampening 

effect  was  described,  for  the  case  of  angularly  misaligned  missiles  during  reentry,  by 
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Figure  25:  Attitude  Dynamics  Coordinate  System 
[Ref.  57] 
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Allen  in  1957.  [Ref.  58]    The  system  stiffness  effect  during  reentry  is  defined  in  the 
quantity 

M[  (94) 

\ 

where 

Mz°    =  the  derivative  of  the  dimensional  aerodynamic  moment  with  respect  to 
the  angle  of  attack 

Iz        =  reentry  body's  moment  of  inertia  relative  to  the  Z-axis 

Based  on  these  facts,  Duzmak  states: 

The  value  of  M2a  increases  by  several  orders  of  magnitude  upon  the  descent 
because  of  the  increase  in  density  and  the  indicated  effect  of  the  variation  in  the 
system's  parameters  is  the  main  factor  determining  the  damping  of  the  oscillation 
...  investigating  the  disturbed  motion  upon  atmospheric  entry  permit  one  conclude 
that  a  descending  entry  body  represents  a  significantly  nonlinear  mechanical 
system.  [Ref.  5 6: p.  7] 

a.    Equations  of  Perturbed  Motion 

The  basic  development  of  the  equations  of  perturbed  motion  of  a  body 
about  its  center  of  mass  during  reentry  is  presented  due  to  the  uniqueness  of  Duzmak' s 
work  [Ref.  56:pp.  12-28].  Equations  (1)  through  (3)  describe  the  system  of  equations 
of  motion  for  this  problem. 

The  characteristic  time  intervals  for  the  reentry  body's  motion  relative 
to  its  center  of  mass  and  the  motion  of  its  center  of  mass  are,   respectively,  defined  as 

'«  ■  ^  (95) 
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t    =  *L  (96) 

V 

where 

0        =  characteristic  rotational  velocity  with  respect  to  the  center  of  mass 
Ar      =  characteristic  variation  in  the  radius  vector  of  the  center  of  mass 
V       =  characteristic  velocity  of  its  center  of  mass 

In  the  development  of  the  equations  of  perturbed  motion,  the  following 

factors  and  assumptions  were  taken  into  account: 


1.  The  interaction  between  the  reentry  body's  motion  with  respect  to  its  center  rf 
mass  and  the  motion  of  the  center  of  mass  is  neglected. 

2.  The  hypothesis  is  satisfied  in  the  upper  atmosphere  portion  of  the  descent 
trajectory  until  the  angles  of  attack  becomes  less  than  one  radian  because  of  the 
atmosphere's  stabilizing  effect. 

3.  The  hypothesis  is  satisfied  in  the  upper  atmosphere  portion  of  the  descent 
trajectory  since  the  body's  kinetic  energy  of  the  center  of  mass  is  many  times 
larger  than  the  body's  kinetic  energy  due  to  its  motion  relative  to  the  center  rf 
mass.  Due  to  this  fact,  the  aerodynamic  moments  begin  to  affect  the  motion 
about  the  center  of  mass  by  stabilizing  the  body  at  significantly  higher  altitudes 
than  the  aerodynamic  forces  acting  on  the  center  of  mass. 

4.  The  center  of  mass  trajectory  parameters  are  known  as  a  function  of  time. 

5.  The  flight  trajectory  is  two-dimensional. 

6.  The  average  rotation  of  the  velocity  vector  of  the  reentry  body's  center  of  mass 
which  is  the  curvature  of  the  average  trajectory  is  neglected.  This  effect  is 
caused  by  the  gravitational  force  and  the  body's  rotational  velocity  about  its 
center  of  mass. 

7.  The  hypothesis  is  clearly  satisfied  for  the  sections  BC  and  DF  as  shown  in 
Figure  26  [Ref.  56:p.l6],  because  the  trajectory  is  nearly  linear  due  to  flight 
velocities  of  several  kilometers  per  second.   In  the  sections  AB  and  FG,  where 
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Figure  26:  Trajectory  Of  The  Reentry  Body 
[Ref.  56]  95 


the  trajectory  inclination  angle,  0,  varies  significantly,  the  hypothesis  is  usually 
satisfied.   The  hypothesis  is  not  satisfied  in  the  exoatmospheric  section  CD. 

8.      The  variation  in  the  velocity  vector  orientation  caused  by  the  effect  of  the 
aerodynamic  lift  force  is  taken  into  account. 


The  reentry  body  coordinate  system,  shown  in  Figure  27  [Ref.  56:p.  23], 
has  both  dynamic  and  aerodynamic  symmetry  and  is  used  for  the  derivation  of  the  three- 
dimensional  perturbed  equations  of  motion: 


1.  The  OXYZ  right-handed  coordinate  system  is  fixed  in  both  inertial  space  and 
relative  to  the  reentry  body. 

2.  The  body's  center  of  mass  is  denoted  by,  0. 

3.  The  OXc  axis  is  the  velocity  vector  of  the  center  of  mass. 

4.  The  angle,  a,  between  the  OX  and  OXc  axis  is  not  the  three-dimensional  angle 
of  attack  or  nutation  angle. 

5.  The  plane  containing  the  OXc,  OX  and  0Y  axis  is  the  angle  of  attack  plane. 

6.  The  motion  of  the  body  in  the  angle  of  attack  plane  is  defined  by  dx/dt  and  is 
directed  along  the  0Z  axis. 

7.  The  rotation  of  the  angle  of  attack  plane  with  respect  to  the  velocity  vector, 
OXc,  is  determined  with  the  help  of  the  precessional  velocity,  X,  which  is 
directed  along  the  OXc  axis. 

8.  The  body's  rotation  with  respect  to  the  angle  of  attack  plane  is  determined  using 
the  intrinsic  rotational  velocity,  /x,  directed  along  the  OX  axis. 

9.  The  angle,  7,  between  the  trajectory  plane  and  the  angle  of  attack  plane  is 
defined  by  d7/dt  =  X. 

10.  The  angles  a  and  7  have  the  following  ranges:  0  <  a  <  180°;  0  <  7  <  00. 


96 


Trajectory   of    the 
\  >»center   of   mass 


Trace   of 
angle   of 


the 


attack's  plane 


6S*r(Q 


Xc 


Trace  of  the 
trajectory' s 
plane 


Figure  27:  Reentry  Body  Coordinate  System 
[Ref.  56] 
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This  coordinate  system  can  also  be  used  to  describe  two-dimensional  motion  about  the 
center  of  mass. 

The  angular  velocity  vector,  co,  of  the  reentry  body  is  the  resultant  of 
three  rotations:  da/dt;  X;  and  /*.  As  previously  mentioned,  the  aerodynamic  lift  force 
affects  the  body's  velocity  vector  orientation.  Lift  acts  in  the  angle  of  attack  plane, 
creating  a  motion  of  the  body's  velocity  vector  that  is  in  the  same  plane  with  the  change 
of  the  angular  velocity  along  the  Z  axis  and  is  defined  as 

Aco ■    -  e  —  (97) 

mV 

where 

L        =  force  of  lift  from  equation  (12) 

m       =  mass  of  the  body 

e        =  t^tv  (small  parameter) 

The  reentry  body's  angular  velocity  vector  X,  Y,  Z  components  are 
defined  as 

w    =  /x  +  X  cos  a  (98) 


wy  =  -X  sin  a  (99) 


u=^£i  (ioo) 

1       dt      mV 

The  reentry  body's  angular  momentum  vector,  H,  XYZ  components  are 
defined  as 
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Hx  =  \wx  (10D 

Hy  =  ^  (102) 

H   =  I  co  (103) 

2  Z      I 

where 

Ixy  z    =  reentry  body's  moments  of  inertia  along  the  X,  Y,  Z  directions 

The  main  moment  force,  M,  acting  on  the  reentry  body  consists  of  a 
restoring  moment,  Mz(t,cx),  and  a  small  damping  moment  that  can  be  projected  in  the  X, 
Y,  Z  directions,  is  given  by 

M(T,a)  =  maAL  (104) 


AO(r,c0  =  m;x^±l  (105) 


Myuy(T,a)  =  m;y^l  (106) 


AO(r,a)  =  mfzim  (107) 


where 


mz  ,  m°x  ,  m°y  ,  /n"z    =  dimensionless  aerodynamic  (108) 

coefficients 

q        =  pV2/2 

T  =   €t 
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By  projecting  dH/dt  =  M  from  equation  (3)  on  the  OXYZ  coordinate 
system  and  taking  equations  (101)  through  (107)  into  account,  the  reentry  body's 
equations  of  perturbed  motion  can  be  written  in  the  form 


da 


(109) 


do)  a 

I  — -  +1  03  w  -I  03  w    -  eM„  y(T,a)o3 

z     jf         x    x    z       z    x    x  j/\')"/wj 


(110) 


d.03 

lz  — -  +  lz  03yxsjx  -  lx  03xwy  =  Mz(t,  a)  +  eM" z(t,  or)  03z 


(111) 


where 


wx       =  Xcosa 


XZy  =     03y 


Equations  (109)  through  (111)  can  be  rewritten  by  substituting  in  equations  (99)  and 
(100)  for  wy,  G7y,  and  wz  and  by  substituting  wx  from  above,  into  the  following  form 


dr 
dt 


+  efx(T,a)r  =  0 


(112) 


dk     (2X  cosa  -r)  da 


dt 


sin  a         dt 


fy(r,a)-——. — 
'  mV  sin  a 


=  0 


(113) 


— r--X2  sin  a  cosa+rX  sina  +  M(t, a) +efz(T, a)  —  =  0 


(114) 
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where 


M(r,a)  = 

Mz(r,a) 

M"jt(t,q;) 

fx(T,a)  = 

I 

(115) 


(116) 


(117) 


/v(r,«)  =  " 


Ky(r,a)  .    L  cos  a 


(118) 


I2  mV  sin  a 


M/z(T,a)     l2 

/z(T,Qt)    = - +  — 

I.  mV 


(119) 


da 


(120) 


Finally,  the  system  of  equations  that  describes  the  two-dimensional 
motion  of  the  reentry  body's  center  of  mass  is  given  by 


dB_ 
dt 


=   € 


CL  cos  y qA  _gT 
mV  V 


cos  6 


1- 


R8r 


(121) 


dV 
dt 


=  c 


C>4 


/w 


-gT  sin0 


(122) 
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^l=eVsmd  (123) 

dt 


£^  =e^Vcos0  (124) 

dt         R 

where 

gT       =  gravitational  acceleration  of  the  Earth 

6        =  angle  between  the  tangent  to  the  trajectory  and  the  local  horizontal 

H       =  altitude  of  the  flight 

R        =  Ro+  H 

L        =  range  of  flight  measured  with  respect  to  the  Earth's  surface 

The  investigation  of  equations  (112)  through  (1 14)  in  the  case  of  known 
solutions  for  equations  (121)  through  (124)  is  the  principal  content  of  Duzmak's  work. 
Specifically,  the  following  major  areas  were  investigated: 


1.  The  distinctive  features  of  the  two-dimensional  motion  that  explains  the 
interaction  of  the  nonlinear  effect  with  the  influence  of  variability  of  parameters. 

2.  The  distinctive  features  of  the  motion  that  results  from  the  transition  through  the 
transonic  velocity  range. 

3.  The  relationship  of  the  angular  momentum  components  with  the  parameters  of 
motion  in  a  increasingly  denser  atmosphere.  This  analysis  clarifies  the  effect 
of  the  shape  of  the  instantaneous  characteristics,  stability  margin  and  damping. 

4.  The  effect  of  the  reentry  body's  motion  parameters  above  the  atmosphere  on  its 
motion  in  the  denser  layers  of  the  atmosphere. 

5.  The  case  of  sinusoidal  dependence  of  the  longitudinal  moment  on  the  angle  of 
attack  for  three-dimensional  motion. 


102 


6.      The  body's  motion  relative  to  its  center  of  mass  for  a  small  angle  of  attack 
reentry  (shallow  angle  reentry  from  a  decaying  orbit). 


b.    Effect  of  Motion  Relative  to  the  Center  of  Mass  on  the  Motion  of  the 
Center  of  Mass 

The  coupled  effect  of  the  perturbed  motion  about  the  center  of  mass  with 
the  trajectory  of  the  center  of  mass  occurs  because  of  the  dependence  of  the  aerodynamic 
coefficients  on  the  angle  of  attack.  Duzmak  neglects  this  interaction  in  the  development 
of  the  asymptotic  solutions.  This  assumption  is  based  on  the  fact  that  the  atmosphere 
will  start  to  influence  the  motion  relative  to  the  center  of  mass  at  higher  altitudes  than 
the  atmosphere  will  start  to  effect  the  motion  of  the  center  of  mass.  Generally,  for 
perturbed  motion  about  the  center  of  mass  in  a  free-molecular  regime,  the  effect  of  the 
lift  force  is  small,  since  it  continually  acts  in  different  directions.  In  the  denser  layer  of 
the  atmosphere,  lift  can  also  provide  some  additional  damping  of  the  angle  of  attack 
oscillation.  [Ref.  56:p.  29] 

The  main  effect  on  the  trajectory  is  exerted  by  adrag,  equation  (6).  For 

perturbed  motion  at  high  altitudes  where  atmospheric  density  is  low  and  CD(a)  varies 

significantly,  a,,,^  has  little  effect  on  the  trajectory.    In  many  cases  at  lower  altitudes, 

perturbed  motion  has  a  weak  effect  on  the  trajectory  due  to  the  fact  that  a  reentry  body 

is: 

...able  up  to  this  instance  to  stabilize  itself  under  the  action  of  aerodynamic 
moments  so  that  the  variation  in  CD  in  the  case  of  perturbed  motion  becomes 
insignificant.  [Ref.  56:p.  29] 


For  the  case  of  hovering  or  similar  motion  where  the  oscillations  do  not 

dampen,  the  oscillations  of  the  angle  of  attack  can  occur  with  significant  amplitude. 

Under  these  conditions,  it  is  essential  to  take  into  account  the  angle  of  attack  oscillation 

effect  on  the  trajectory  of  the  center  of  mass.  [Ref.  56:p.  317] 

The  determination  of  the  range  component  of  the  scattering  landing  points  is 
important.  This  scattering  is  caused  by  the  oscillations,  the  angle  of  attack  and  the 
variation  of  the  coefficient  of  drag  associated  with  the  oscillations.  This  quantity 
cannot  in  general  be  determined  without  taking  into  account  the  interaction  of 
motion  relative  to  the  center  of  mass  and  the  motion  of  the  body  center  of  mass. 
[Ref.  56:pp.  317-318] 

In  the  derivation  of  a  simple  approximate  method  to  solve  this  problem, 

the  descent  trajectory  was  divided  into  two  sections:   H    >    Ht,  (trajectory  in  the 

atmosphere's  tenuous  layers)  and  H  <  Hb  (trajectory  where  the  interaction  between  the 

motion  with  respect  to  the  center  of  mass  and  the  center  of  mass  trajectory  is  taken  into 

account).    Hb  is  the  limiting  altitude  that  is  calculated  to  be  70-80  km..   In  the  region, 

H  <  Hb,  the  refined  asymptotic  method  breaks  down.    An  approximate  solution  is 

derived  by  formulating  averaged  equations  for  the  asymptotic  solutions  that  describe  the 

change  of  the  slow  varying  components  of  : 

1.  The  maximum  and   minimum  values  of  the  angle  of  attack  during  each 
oscillation 

2.  The  instantaneous  oscillation  period 

3.  The  averaged  components  of  the  center  of  mass  trajectory 

The  derivation  of  this  method  is  mathematically  rigorous.  [Ref.  56:pp.  319-336] 
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The  direct  solution  to  this  problem  can  be  obtained  by  numerically 
integrating  the  complete  equations  of  motion  given  by  [Ref.  56:p.  45] 


d2a+MJ.(?cl^1...ta)+effic1jc2,....,ct)-£  =  0 


(125) 


— -+eXi(xljc2,...ia)  =  0 
dt 

a  =  u,...,6) 


(126) 


where 


xm         (G-r  cosa)(r-G  cosa)     x.,      s 
MT  =  -i Jl L  +M(r,a) 

sin'1  a 


(127) 


Xj  =fy(T,a)G  +  ([fx(T,a)-f(T,a)]  cos  a  -L{t,o)  sina)r 


X2  -fxij,a)r 


(128) 
(129) 


CJa)qA 


m 


(130) 


X.  = 


gT  cos  6 


1- 


r8t 


(131) 


X5  =  -V  sin0 


(132) 


X,  =  -_2Vcos0 
6         R 


(133) 
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^       Ocosa-n   sine* 

G  =  — y- =  A  siir  a  +  r  cost*  =  constant  (134) 

However,  in  1970,  this  numerical  integration  required  an  excessive  amount  of  machine 
time  or  computer  time,  since  the  solution  contains  high  frequency  components  along  the 
trajectory. 

c.     Follow-on  Investigations  of  Uncontrolled  Reentry  Body  Motion 

References  [59]  and  [60]  are  extensions  of  Duzmak's  investigation.  Each 
paper  examines  a  certain  aspect  of  the  uncontrolled  three-dimensional  motion  of  a  reentry 
body  relative  to  its  center  of  mass. 

Reference  [59]  investigates  the  three-dimensional  uncontrolled  motion  of 
a  spherical  body  relative  to  its  center  of  mass  with  an  arbitrary  angle  of  attack.  The 
mean  differential  equations  are  derived  for  rolling  motion.  These  equations  are  solved 
in  implicit  form  for  any  angle  of  attack  by  using  the  Van  der  Pol  method  of  integration. 
Reference  [60]  investigates  the  three-dimensional  uncontrolled  motion 
relative  to  the  center  of  mass  of  a  reentry  body  with  a  small  geometric  and  dynamic 
asymmetry.  An  approximate  analytical  solution  for  the  equations  for  unperturbed  three- 
dimensional  and  the  averaged  equations  for  perturbed  three-dimensional  motion  are 
derived.  By  imposing  no  limits  on  the  initial  angular  velocity  and  the  three-dimensional 
angle  of  attack  or  nutation  angle,  the  averaged  equations  of  motion  computational  time 
is  reduced  by  a  factor  of  approximately  three  as  compared  to  the  numerical  integration 
of  the  complete  equations  of  motion. 


106 


B.      SIX-DEGREE-OF-FREEDOM  SIMULATIONS 

The  foundation  for  six-degree-of-freedom  motion  has  been  discussed  in  the  previous 
section.  Several  other  investigators  have  examined  the  effect  of  motion  relative  to  the 
center  of  mass  on  the  trajectory  of  the  center  of  mass  [Refs.  61-64].  The  equations  of 
orbital  decay,  (121)  through  (124),  when  coupled  with  the  equations  of  attitude  motion, 
(112)  through  (114),  will  completely  describe  the  motion  of  a  vehicle  during  the  final 
stages  of  life.  Equations  (125)  and  (126)  are  the  coupled  equations  of  six-degree-of- 
freedom  motion  written  as  a  function  of  angle  of  attack. 

Coupling  these  equations  means:  the  solution  of  one  set  of  equations  determines  the 
magnitude  and  form  of  the  forces  or  moments  in  the  other  set  [Ref.  62 :p.  13].  For 
example,  if  the  area  changes  due  to  either  a  loss  in  mass  or  a  change  in  attitude,  then  the 
ballistic  coefficient  (W/CDA)  and  the  coefficient  of  drag  will  be  affected. 

References  [61],  [62]  and  [63]  investigate  Skylab's  attitude  and  trajectory  motion. 
Specifically,  references  [62]  and  [63]  derive  a  variation  of  the  coupled  equations 
presented  in  the  previous  section.  The  authors  performed  a  six-degree-of-freedom 
trajectory  simulation  using  the  coupled  equations  of  motion.  Figure  28  shows  the  flow 
chart  used  in  their  simulation  [Ref.  62:p.  14].  Numerical  integration  of  the  six 
differential  equations  yielded  altitude  and  orbital  elements  as  a  function  of  time. 
Computer  run  results  for  these  simulations  were  extremely  long  and  therefore  the 
decision  was  made  to  artificially  increase  the  magnitude  of  the  aerodynamic  damping. 
Subsequently,  the  results  do  not  simulate  the  actual  dynamical  behavior,  but  they  show 
the  "possible"  types  of  dynamical  behavior  for  Skylab.  [Ref.  62:p.  15] 
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Reference  [64]  is  a  reentry  analysis  of  a  proposed  radioisotopic  thermoelectric 
generator  (RTG)  that  was  connected  to  a  generic  user  satellite.  A  3-degree-of-freedom 
reentry  trajectory  simulation  was  conducted  that  utilized  aerodynamic,  material  property 
and  heating  characteristics. 

Newtonian  and  free  molecular  drag  coefficients  were  calculated  using  the  Mark  IV 
Supersonic-Hypersonic  Arbitrary-Body  program  for  the  generic  satellite  and  each  of  its 
subelements.  Figure  29  shows  the  three-dimensional  shapes  of  the  generic  spacecraft  and 
one  subelement  (reactor  with  radiator)  used  by  the  program  [Ref.  64:p.  6].  Aerodynamic 
blockage  was  neglected  and  a  zero  angle  of  attack  was  assumed  in  the  calculations.  The 
heat  transfer  methodology  used  in  the  simulation  to  predict  the  satellite's  breakup  during 
reentry,  implemented  a  calibrated  heat  transfer  model  to  closer  simulate  actual  conditions 
[Ref.  64:pp.  7-12].  This  will  be  discussed  further  in  the  section  on  breakup  later  in  the 
chapter.  The  trajectory  simulation  was  designed  to  change  the  mass  properties  and 
aerodynamic  coefficients  as  the  shape  changed  due  to  an  "assumed  predetermined 
breakup"  scenario.  This  assumption  was  based  on  the  basic  geometric  components  or 
subelements  of  the  satellite  and  their  probable  separation  sequence  during  breakup.  For 
example,  the  boom  separated  from  the  main  part  of  the  spacecraft  and  the  reactor  failed 
first,  this  was  followed  by  the  heat  radiator  cone,  and  then  the  other  components  in 
sequence.  [Ref.  64: p.  14] 

Trajectory  simulations  were  run  to  assess  the  breakup  altitudes,  due  to  reentry 
heating,  and  minimum  footprint  lengths,  due  to  fuel  pin  release  at  various  altitudes, 
which  were  independent  of  the  heating  effects. 


109 


GENERIC  SPACECRAFT 


REACTOR  WITH  RADIATOR 


Figure  29:  Simulated  Generic  Spacecraft  And  Reactor  Sub-Element 
[Ref.  64] 
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C.      RAREFIED  GAS  DYNAMICS 

Numerous  authors  have  investigated  reentry  vehicle  attitude,  heating  rates  and 
aerodynamic  coefficients  in  relation  to  the  atmospheric  flow  regimes  [Refs.  22,  56,  65- 
75].  As  stated  in  Chapter  II,  the  changes  in  flow  regime  and  corresponding  changes  in 
critical  parameters  of  the  reentry  trajectory  or  heating  equations,  are  poorly  understood 
in  the  transition  from  free-molecular  to  hypersonic  continuum  flow.  For  engineering 
applications  the  quantities  of  lift,  drag  and  heat  transfer  are  usually  estimated  by  a 
"judicious  faring"  between  regimes   [Ref.  22 :p.  203]. 

One  solution  to  this  problem  is  presented  by  Koppenwallner  and  Johannsmeier, 
reference  [76].  This  solution  is  a  technique  of  "bridging"  between  the  free  molecular  and 
continuum  flow,  based  on  Newtonian  and  free  molecular  theory.  This  technique  allows 
the  prediction  of  lift  and  drag  during  the  hypersonic  entry  phase.  [Ref.  76: p.  461] 

Three  hypersonic  flow  regimes,  with  approximate  boundaries,  are  described  as 
follows:  [Ref.  76:p.  461] 

1.  Free  molecular  flow  Kn  >  5 

2.  Rarefied  transitional  flow  5  >  Kn  >  0.001 

3.  Hypersonic  continuum  flow  Kn  <  0.001 

These  boundaries  are  actually  dependent  upon  reentry  vehicle  shape  and  on  the 
aerodynamic  property  considered.  Figure  30  shows  these  flow  regimes  in  an  altitude- 
velocity  graph  [Ref.  76:p.  461]. 
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Figure  30:  Flow  Regimes  -  Altitude  vs  Velocity 
[Ref.  76] 
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The  analysis  of  this  technique  is  limited  to  flow  conditions  where  the  hypersonic 
Mach  independence  principle  is  applicable;  therefore,  blunt  shaped  reentry  vehicles  will 
have  a  lower  limit  of  Mach  5.  Figure  31  shows  the  typical  aerodynamic  data  variation, 
for  a  blunt  body,  through  the  three  flow  regimes  of  interest  [Ref.  76:p.  461]. 

The  typical  drag  coefficient  behaves  such  that  in  free  molecular  flow,  the  drag 
coefficient  is  independent  of  Knudsen  number  and  the  value  of  CD  =  2.  The  drag 
coefficient  decreases  throughout  the  transitional  flow  regime  and  again  reaches  a  constant 
value  in  the  continuum  flow  regime.  [Ref.  76:p.  461] 

The  typical  lift  coefficient  behaves  such  that  in  the  free  molecular  flow,  the  lift  is 
very  small  or  negligible.  The  lift  coefficient  increases  throughout  the  transitional  flow 
regime  until  it  reaches  a  hypersonic  continuum  flow  value.  [Ref.  76:p.  461] 

The  heat  transfer  Stanton  number  (St)  is  very  close  to  a  value  of  one  in  the  free 
molecular  flow  regime.  In  the  transitional  flow,  the  heat  transfer  efficiency  decreases 
and  the  Stanton  number  decreases.  In  the  continuum  flow,  the  Stanton  number  decreases 
continuously  as  a  function  of  Reynolds  number  in  the  manner  St—  l/v^Re.  [Ref.  76:p. 
461] 

Table  III  describes  the  bridging  dependencies  as  modeled  above.  The  derived 
transitional  functions  D2  and  L2  are  functions  of  angle  of  attack,  vehicle  shape  and 
Knudsen  number.  These  transitional  functions  must  bridge  the  free  molecular  and 
continuum  flow  regimes  and  must  degenerate,  in  the  two  limits,  to  the  free  and 
continuum  flow  values  as  shown  in  Table  IV.  [Ref.  76] 
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Figure  31:  Aerodynamic  Coefficients  In  The  Flow  Regimes 
[Ref.  76] 
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Table  III: 
[Ref.  76] 


BRIDGING  DEPENDENCIES 


FLOW  REGIME 

DRAG 

LIFT 

Free  molecular 

CD  =  D,  (a,  shape) 

CL  =  0 

Transitional 

CD  =  D2  (a,shape,  Kn) 

CL  =  L2  (a,shape,Kn) 

Hypersonic 
continuum 

CD  =  D3  (a,  shape) 

CL  =  L3  (a,  shape) 

Table  IV: 
[Ref.  76] 


TRANSITIONAL  BRIDGING  FUNCTIONS 


DEGENERATE  TO 
FREE  MOLECULAR 

TRANSITIONAL 
BRIDGING  FUNCTION 

DEGENERATE  TO 
HYPERSONIC  CONT. 

D,  =  {a,  shape) 

D2  =  (a,  shape,  Kn) 

D3  =  {a,  shape) 

L,  =  0 

L2  =  (o,  shape,  Kn) 

L3  =  (a,  shape) 

Local  flow  independence  exists  in  free  molecular  as  well  as  hypersonic  Newtonian 
flow.  This  means  that  unless  "shadowing"  exists,  surface  and  shape  elements  will  not 
influence  each  other.  The  shape  elements  of  this  technique  are  cones,  spherical  caps  and 
cylindrical  shells.  These  shape  elements  allow  the  composition  of  capsules,  blunted 
cones  and  cone-cylinders.  Figure  32  shows  the  shape  elements  and  several  composed 
bodies  [Ref.  76:p.  462]. 
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Figure  32:  Shape  Elements  And  Composed  Bodies 
[Ref.  76] 
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1.       Newtonian  Aerodynamics  for  Hypersonic  Continuum  Flow 

For  the  case  of  blunt  bodies,  it  can  be  assumed  that  the  main  contribution  to 
drag  is  pressure;  therefore,  the  influence  of  viscous  effects  on  the  aerodynamic 
coefficients  are  neglected.  [Ref.  76:p.  462]  Using  the  Newtonian  pressure  law,  the 
coefficients  of  drag  and  lift  may  be  described  in  differential  equations.  The  basic 
Newtonian  relationships  are  as  shown  in  Table  V.  [Ref.  76:p.   462] 


Table  V: 

[Ref.  76] 


NEWTONIAN  LOCAL  PRESSURE  LAW 


Wetted  surface 

6    <  90° 

p  /  q06=  kN  cos20 

Newtonian  wake 

6    >  90° 

p  /  q«,  =  CDo 

where 


P 

q*.     = 

^       = 

e 

K  = 


pressure 

dynamic  pressure  in  free  stream 

Newton  factor 

flow  inclination  against  surface  normal 

ratio  of  specific  heats 

coefficient  of  drag  as  a  function  of  angle  of  attack 


and 

ktf      =2  (simple  Newton) 

kw      =  (k  +  3)/(k  +  1)  (modified  Newton) 

In  the  Pike  formulation,  the  Newtonian  differential  equation  is  given  by  [Ref.  71  :p.  462] 

CD(a)  *  CD(a)cot(a)  *  12CD(a)  =  6kN  ^L  (135) 

AR 

3CL(a)  =  CDm  (136) 

where 

Ap      =  Flow  projected  area  of  body 

AR      =  Reference  area  (7rD2/4) 

In  order  to  determine  the  drag  coefficient  and  solve  the  differential  equations, 
the  following  approach  is  taken:  CD  (a=0°,  shape).  This  implies  that  the  zero  angle 
of  attack  drag  is  proportional  to  the  Newton  factor  and  to  a  shape  dependent  factor,  Cs 
Ref.  76:p.  462] 

CD(a=0°)  =  kN  Cs  (137) 

Upon  integration  of  the  Newtonian  pressure  distribution,  the  following  values 
for  Cs  are  determined  for  various  bodies  in  Table  VI.  [Ref.  76:p.  462] 
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Table  VI: 

[Ret.  76] 


SHAPE  FACTOR  VALUES 


Disk 

cs=  1.0 

Sphere 

Cs=  0.5 

Cylinder 

Cs=  0.67    (cross  flow) 

Cylindrical  shell 

Cs=  0          (parallel  flow) 

Sharp  cone 

Cs=  sin20 

Spherical  cap 

Cs=  1  -  1/8  •  (d/rN)2 

Figure  33  shows  the  shape  factor  as  a  function  of  the  ratio  (d/rN),  diameter 
over  nose  radius  [Ref.  76:p.  466]. 

When  the  shape  factor  or  drag  at  zero  angle  of  attack  are  known,  the 
aerodynamic  coefficients  of  the  class  are  completely  fixed.  The  following  universal 
solutions  for  drag  and  lift  are  obtained  from  the  Newtonian  differential  equations  [Ref. 
76:p.  462] 
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Figure  33:  Shape  Factor  As  A  Function  Of  (d/rN) 
[Ref.  76] 
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CD  =  JL  (2CS  -  (5CS  -  3)  sin2(c*))  cost*  (138) 

CL  -  4*  (  2  (1  -  2CS)  -  (5CS  -  3)  sin>(a))  sine*  (139) 

These  solutions  are  valid  for  a  <  a,,^  ,  since  at  a^  the  wetted  area  decreases  due  to 
Newtonian  shadowing.  Figure  34  shows  the  universal  Newtonian  lift  and  drag  functions 
for  various  angles  of  attack  [Ref.  76: p.  466]. 

The  shape  factor  at  zero  angle  of  attack  serves  as  the  critical  parameter  in  this 
method.  The  solution  is  valid  only  under  the  condition  that  the  wetted  surface  area 
remains  constant  for  any  angle  of  attack.  [Ref.  76:p.  462] 

As  the  geometric  body  bluntness  increases,  Cs  increases  and  Newtonian 
shadowing  is  shifted  to  higher  angles  of  attack.  Table  VII  shows  that  depending  upon 
the  body  shape,  the  lift  slope  at  zero  angle  of  attack  may  be  either  positive  or  negative. 
[Ref.  76:p.  462] 
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Figure  34:  Newtonian  Lift  And  Drag  Functions 
[Ref.  76] 
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Table  VII: 

[Ref.  76] 


BODY  SHAPE  AND  LIFT  SLOPE 


Cs 

Body  shape 

Lift  slope  @  a  =  0 

<  0.5 

slender 

positive 

>  0.5 

blunt 

negative 

This  demonstrates  that  blunt  reentry  bodies  will  experience  a  negative  lift  for  a  positive 
defined  angle  of  attack.  [Ref.  76:p.  462] 

2.       Free  Molecular  Flow  Model 

The  technique  of  reference  [76]  also  develops  analytical  formulas  for  free 
molecular  flow.   The  general  considerations  include:  [Ref.  76:p.  463] 

1.  Accommodation  coefficient,  a 

2.  Finite  molecular  speed  ratio,  S 

3.  Wall  temperature,  Tw 

These  factors  alone  are  insufficient  thus  the  following  simplifying  assumptions  are  made: 
[Ref.  76:p.   463] 


1.  Body  shape 

2.  Aerodynamic  cold  wall 

c 

3.  Diffuse  molecular  reflection 


Hypersonic  blunt,  Ma  (d/1)  >  2 
Tw/T,,    <<  1 
a  =  1 
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From  these  constraints  the  local  surface  pressure  and  shear  stress  may  be 

defined  as  shown  in  Table  VIII.  The  aerodynamic  force  coefficients  are  also  derived  as 

shown  in  Table  IX. 

Table  VIII:  LOCAL  SURFACE  PRESSURE  AND  SHEAR  STRESS 
[Ref.  76] 


Wetted  surface 

G    <  90° 

p/q0O=  2  cos20 

Unwetted  surface 

6    <  90° 

r/q00=  sin0«cos0 

Table  IX: 

[Ref.  76] 


AERODYNAMIC  FORCE  COEFFICIENTS 


Drag  coefficient 

CD=  2«AP(flr)  /  AP(0)  =  cos  a 

Lift  coefficient 

CL=  0 

These  assumptions  imply  that  no  lift  is  produced  and  drag  is  proportional  to  the  flow 
projected  area,  AP(a),  in  the  free  molecular  flow  regime.  [Ref.  76:p.  463] 

In  the  more  general  free  molecular  flow  case,  where  Tw/T^  ,  an  and  a,  are 
considered,  an  axisy  metric  blunt  body  formulation  for  the  force  coefficients  may  be  given 
by 
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c„  = 


zacosa  +  a    - — 

1  n      s 


W 


c»*i 


2      si 


cos2a 


(2-a„-a,U[3(C52+l)cosa+(5C52-3)cos3a 


(140) 


CI                           V   «                   w 
,    -   o_- —     


yx 

T 

> 

OB 

c  -1 

S1    2 


sin2a  + 


(2-an-a;)i[-(C52+l)sina-(5Cs2-3)sin3a] 


(141) 


where 

To,      =  temperature,  free  stream 

CS1     =  shape  coefficient,  in  front  of  normal  shock  wave 

Cs2     =  shape  coefficient,  behind  normal  shock  wave 

The  physical  significance  of  several  of  the  terms  in  equations  (140)  and  (141) 
are:  [Ref.  76:p.  463] 


1.  The  first  term  in  CD  gives,  for  diffuse  reflection,  a=\,  the  contribution  of  the 
incident  flux  to  the  aerodynamic  coefficients. 

2.  The  Tw/Tw  term  states  the  influence  of  the  reentry  vehicle  surface  wall 
temperature  on  drag  and  lift. 

3.  The  (2-cvaJ  term  vanishes  for  diffuse  molecular  reflection. 


In  simple  free  molecular  flow  the  equations  will  degenerate  appropriately  as 
aD  =  1,  <Tt=  1  and  Tw/T^  =0.  In  the  Newtonian  formulation  the  equations  will  degenerate 
as  an=l,  at=0  and  Tw/To^O.  In  this  last  case  Csl  vanishes  and  Cs2  is  the  Newtonian 
shape  coefficient,  Cs.    [Ref.  76:p.  463] 
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3.       Bridging  Free  Molecular  and  Continuum  Flow 

Several  methods  exist  which  attempt  to  bridge  the  gap  between  the  free 
molecular  and  continuum  flow  regimes.  [Refs.  77-79].  In  the  former  USSR  and  DLR, 
local  bridging  is  accomplished  with  a  finite  surface  element  method,  [Refs.  77-78].  In 
the  U.S.,  reference  [79],  bridging  has  been  accomplished  through  an  integral  coefficients 
method.    [Ref.  76:p.  463] 

The  method  of  bridging  by  shape  element  description  is  developed  in 
reference  [76]  and  is  presented  in  the  following  section.  The  methods  of  references  [78] 
and  [79]  may  be  used  to  derive  analytical  formulas  for  trajectory  calculation.  However, 
the  basic  bridging  relations  must  be  derived  through  experimentation.  [Ref.  76:p.  463] 

a.     Shape  Element  Bridging  Method. 

Experimental  data,  presented  in  Figures  35  and  36  show  the  drag 
coefficient  changes  of  a  sphere  and  disk  respectively.  The  data  covers  the  the  entire 
transitional  flow  regime  [Ref.  76]. 

For  a  sphere,  the  Reynolds  number,  behind  the  normal  shock  (Re^),  can 
be  related  to  the  Knudsen  number,  in  the  free  stream  (Kni),  as  follows  [Ref.  76:p.  463] 


'2 


2* 

'     1     ' 

Re,  =  1.26    JL    J-  (142) 


From  this  experimental  data,  an  approximate  formula  for  CD,  as  a 
function  of  Reynolds  number  and  shape,  can  be  derived  by  using  the  "reduced" 
aerodynamic  coefficients.  [Ref.  76:p.  463] 
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Figure  35:  Sphere  Drag  Coefficient  In  Rarefied  Flow 
[Ref.  76] 


127 


4iea 
nn  3034 

2025  1Q135 


300-  5C0 
3CO-5CO 

6C0 

300 


65         f'e*  »»< 

1C9  12      IV  con  nor 
J12-D1     \:%'3*  nox 
8      20  i  *"'  .  i« 


QCIv    tXJt 

gatv  bol 
StrQ.ngauQ* 
»i  maqn  bat 


OOt 


Ql 


10 


10 


KO 


BOO 


Rfj 


Figure  36:  Disk  Drag  Coefficient  In  Rarefied  Flow 
[Ref.  76] 
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_  Cd(^-Cdc 
C        -C 

Z725  m        L/D(a=const.,  Re2)  (144) 

LID  (a  ■  const. ,  Continuum) 

where 

Cdc    =  drag  coefficient,  continuum 

Qjfm  =  drag  coefficient,  free  molecular 

CD      =  reduced  drag  coefficient 

L/D    =  reduced  lift-drag  ratio 

Analysis  of  the  reduced  coefficients  at  zero  angle  of  attack  shows  the 
following:  [Ref.  76:p.  463] 

1.  Slope  of  the  reduced  drag  coefficient  is  shape  independent. 

2.  Continuum  and  free  molecular  flow  boundaries  are  shape  dependent. 

3.  Continuum  and  free  molecular  drag  coefficients  are  shape  dependent. 

Therefore  the  bridging  function  is  derived  from  the  reduced  coefficients  as  follows  [Ref. 
76:p.  464] 

1  Re 

Cn  =  — —  (CDFU  -  Cnr)  In  _ i£  +  Cnc  (145) 

D      5.205     DFM      DC)      Re2       DC 

And  from  this,  the  drag  coefficient  may  be  extracted  in  one  of  two  forms 
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C°  ""  J^05(CDFM'Cp^lnTc+CDC 


(146) 


cD- 


m  CD(Re2)-CD  m      i    inRe2c 


CDFM~CDC  2'26         Re2 


(147) 


The  shape  dependent  boundaries  are  defined  as: 

1.  Re^c  >  Re2  >  RejpM 

2.  Knc  <  Kn  >  Kn™ 

and  the  experimentally  derived  values  are  summarized  in  Table  X  below.  [Ref.  76:p. 
464] 


Table  X: 
[Ref.  76] 


SHAPE  DEPENDENT  BOUNDARIES 


Body 
Shape 

d/RN 

Continuum 
Re2C 

Flow 
Knc 

Free  Molecular  Flow 
Re2FM                KnFM 

Disk 

0 

28.5 

0.117 

0.115                   21.5 

Sphere 

2 

89.3 

0.037 

0.492                      6.8 

In  summary,  the  technique  of  reference  [76]  allows  the  following 
conclusions  to  be  drawn:  [Ref.  76:p.  464] 


1.  Newtonian  theory  based  on  Pike's  method  for  the  hypersonic  continuum  flow 
is  useful,  as  it  shows  the  shape  dependent  and  shape  independent  aerodynamic 
coefficients. 
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2.  The  bridging  function  is  based  on  experimentally  derived  results  for  the 
transition  flow. 

3.  Lift  and  drag  bridging  do  not  follow  the  same  laws. 


b.    Local  Bridging  Method 

Reference  [78]  develops  a  technique  of  "local"  bridging.  The  differences 
between  the  shape  element  method,  referred  to  as  a  "global"  method,  and  local  bridging 
are:  [Ref.  78:p.  469] 


1.  Global  bridging  performs  across  the  spectrum  of  aerodynamic  coefficients  for 
the  complete  body. 

2.  Global  bridging  is  usually  tailored  to  a  specific  class  of  shapes. 

3.  Global  bridging  techniques  require  new  experimentally  derived  "fitting 
constants"  for  each  new  shape  of  interest. 

4.  Local  bridging  is  a  method  of  bridging  the  local  pressure  and  shear  stress 
coefficients.  The  local  distribution  is  then  integrated  over  the  body  surface 
which  yields  the  global  force  and  moment  coefficients. 


The  common  principle  of  all  bridging  techniques  is  the  manner  in  which  they  model  the 
transition  function.  The  known  free  molecular  and  continuum  flow  limiting  values  for 
a  specific  aerodynamic  coefficient  (lift,  drag  or  heat  transfer)  are  weighted  and  applied 
to  the  bridging  function.   The  general  case  is  as  follows 

C(X)  =  CFMfiX)  ♦  Cc  (1  -AX))  (148) 


where 

C        =  aerodynamic  coefficient  considered 

X       =  rarefaction-dependent  flow  parameter 

The  potential  improvements  of  local  bridging  over  global  bridging  are: 
[Ref.  78:p.  470] 

1.  Ease  of  adaptation  to  different  shapes  without  changing  internal  constants. 

2.  More  reliable  calculation  of  moments.  These  are  very  sensitive  to  rarefaction 
because  of  the  varying  contributions  of  pressure  and  shear  in  the  different  flow 
regimes. 

3.  Simplified  method  for  determining  reference  quantities  for  local  coefficients, 
such  as  reference  length  based  on  a  local  coordinate. 

4.  Ability  to  account  for  non-uniform  flow  field  conditions. 

The  conclusion  of  reference  [78],  after  examining  several  different 
bridging  methods  and  various  shapes,  is  that  the  current  state  of  experience  in  applying 
the  local  bridging  methods  allows  no  definitive  answer  as  to  whether  or  not  they  can 
serve  as  an  effective  analysis  tool  in  transition  flow  analysis. 

4.       Gas-Surface/Gas-Gas  Interactions 

Gas-surface  interactions  are  significant  in  the  understanding  of  reentry 
dynamics  and  aerothermodynamics.  The  primary  cause  of  concern  is  that  at  orbital 
altitudes,  the  highly  rarefied  flowfield  is  dominated  by  gas-surface  interactions  that  occur 
at  average  velocities  corresponding  to  that  of  the  orbiting  vehicle.  [Ref.  80:p.  1]  Under 
these  reentry  conditions,  gas-gas  interactions  become  important  and  gas  molecules  that 
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reach  the  vehicle  surface  tend  to  have  lost  some  of  their  initial  translational  energy.  This 
loss  of  translational  molecular  energy  is  due  to  conversion  into  other  forms  such  as  heat, 
internal  energy,  etc.  Figure  37  shows  molecular  velocity  distributions  in  rarefied  and 
transitional  flow  about  a  reentering  sphere  [Ref.  80].  The  collision  mechanisms  by 
which  the  translational  energy  is  converted  include: 

1.  Chemical  reactions 

2.  Ionization 

3.  Dissociation 

Since  the  nature  of  the  gas-surface  interaction  is  known  to  be  dependent  upon  the  velocity 
and  energy  of  the  incident  molecules,  it  becomes  necessary  to  know  the  state  of  the  gas 
molecules  reaching  the  surface.  [Ref.  80:p.  1] 

Wilmoth,  et  al,  [80]  uses  the  technique  of  Direct  Simulation  Monte  Carlo 
(DSMC)  as  developed  by  Bird,  [Ref.  81]  where  the  molecular  velocity  and  energy 
distributions  of  the  gas  molecules  are  a  direct  result  of  the  simulation  process.  Bird  uses 
a  simple  engineering  model  of  the  gas-surface  interactions,  which  accounts  for  diffuse 
and  specular  reflection  along  with  other  phenomena.  This  model  can  accommodate 
processes  such  as  catalytic  reactions  and  molecular  recombination.  A  limitation  of  this 
method  is  that  parametric  studies  must  often  be  performed  in  order  to  place  bounds  on 
the  predicted  quantities  of  interest.  [Ref.  80:p.  1]  This  also  implies  that  the  analyst  must 
judiciously  apply  the  model  based  on  experience  and  a  limited  base  of  experimental  data. 
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Figure  37:  Molecular  Velocity  Distribution 
[Ref.  80] 
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The  lack  of  experimental  data  is  the  reason  why  more  detailed  gas-surface 
interaction  models  have  not  been  developed.  This  results  mainly  for  two  reasons:  [Ref. 
80:p.  1] 


1.  It  is  very  difficult  to  simulate  gas-surface  interactions  at  orbital  or  entry 
velocities  in  a  laboratory. 

2.  It  is  difficult  to  characterize  the  surfaces  used  in  laboratory  experiments  with 
sufficient  generality  that  the  results  may  have  application  in  an  engineering 
context. 


Because  of  the  tendency  of  gas  molecules  to  decelerate  after  gas-gas  collisions,  which 
occur  before  reaching  the  vehicle  surface,  it  is  important  to  quantify  the  actual  velocities 
and  energies  encountered  in  such  cases. 

Reference  [80]  studies  the  "typical"  gas-surface  interactions  for  transitional 
flow  at  entry  velocities.  The  reentry  vehicle  was  a  1.6  m  diameter  sphere  at  a  free 
stream  velocity  of  7.5  km/sec  over  an  altitude  range  of  90  to  130  km.  [Ref.  80:p.  2] 
The  study  was  conducted  using  the  DSMC  method  of  Bird. 

The  flow  conditions  of  the  DSMC  simulation  are  given  in  Table  XI  [Ref. 
80:p.  6].  The  atmosphere  is  modeled  by  Jacchia,  1977,  with  an  exospheric  temperature 
of  1200K.  The  surface  temperature  of  the  sphere  was  assumed  constant  at  350K.  The 
gas-surface  interaction  was  assumed  to  be  diffuse  with  full  thermal  accommodation,  and 
the  surface  was  non-catalytic.  Five  atmospheric  molecular  species  were  modeled  02,  N2, 
O,  N  and  NO,  with  23  reaction  possibilities.  [Ref.  80:p.  3] 
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Table  XI: 
[Ref.  80] 


FREESTREAM  CONDITIONS 


Altitude 
km 

kg/m3 

v_ 
km/s 

T. 
K 

Mole 
Fraction 

Mole 
Fraction 
N, 

Mole 
Fraction 
0 

M 
g/mol 

-*- 
m 

90 

3.43x10" 

7.5 

188 

0.209 

0.788 

0.004 

28.80 

0.017 

100 

5.66x10' 

7.5 

194 

0.177 

0.784 

0.040 

28.24 

0.100 

110 

9. 67x10* 

7.5 

247 

0.123 

0.770 

0.106 

27.22 

0.599 

120 

2.27x10* 

7.5 

368 

0.085 

0.733 

0.183 

26.14 

2.681 

130 

8.23x10* 

7.5 

500 

0.071 

0.691 

0.238 

25.43 

7.724 

Wilmoth's  results  show  the  following:  [Ref.  80:pp.  3-5] 


1.  At  130  km,  a  density  rise  of  nearly  22  times  the  freestream  value  occurs  over 
a  distance  of  —  3  m.  However,  based  on  analysis  of  the  data  it  is  determined 
that  very  few  collisions  are  occurring. 

2.  At  90  km,  the  density  increases  to  well  over  100  times  the  freestream  value  over 
a  distance  of  —0.1  m.  In  this  altitude  range  collisions  become  very  significant 
and  there  is  considerable  chemical  activity. 

3.  There  are  significant  variations  in  the  velocity  and  energy  of  molecules  reaching 
the  surface  over  an  altitude  range  of  130  to  90  km.  Figure  38  shows  the 
average  translational  energy  per  particle  striking  the  reentry  vehicle  surface 
[Ref.  80]. 


D.      SURFACE  ROUGHNESS  EFFECTS 

Another  important  consideration  in  the  reentry  phase  is  the  "boundary-layer 
transition."  It  is  this  phenomenon  which  is  responsible  for  heat  transfer  to/from  the 
reentry  body  due  to  atmospheric  contact  with  the  body.  In  the  continuum  flow  regime 
viscous  effects  are  essentially  restricted  to  a  small  layer  called  the  boundary  layer.   It  is 
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Figure  38:  Average  Translational  Energy 
[Ref.  80] 
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within  this  boundary  layer  that  the  details  of  "fluid"  motion  determine  the  levels  of  skin 
friction  and  heat  transfer  from  the  flow.  [Ref.  22:p.  211] 

A  parameter  commonly  used  for  characterizing  the  boundary  layer  is  the  Reynolds 
number.  [Refs.  82-83]  The  Reynolds  number  is  defined  as  the  ratio  of  inertial  forces  to 
the  viscous  forces,  as  given  by  [Ref.  22 :p.  212] 

inertia  forces    m  A  (momentum) I  (unit  time)    :R  (\49) 

viscous  forces       (shear  stress)x(unit  area) 

In  order  for  a  fluid  (the  atmosphere  is  modeled  as  a  fluid  in  the  continuum  flow) 
to  support  a  shear  stress  there  must  be  relative  motion  between  adjacent  layers.  This 
implies  that  there  is  a  velocity  differential  within  the  flow  layers.  The  following  terms 
are  defined:  [Ref.  22:p.  213] 

t        =  shear  stress  =  ^(dV/dy) 

/x        =  dynamic  viscosity 

y        =  coordinate  normal  to  direction  of  motion 
The  definition  of  Reynolds  number  may  now  be  given  as  [Ref.  22:p.  213] 

Re  m    d(mV)ldt    m  [PL'VI(LIV)\  m     VL  m  VL  (150) 

fi(dV/dy)A  fi(V/L)L2  fi  v 

where 

v        =  kinematic  viscosity  =  (ply.) 

Now  the  Reynolds  number  may  be  used  to  characterize  the  boundary  layer  flow 
conditions.  When  the  viscous  forces  are  large  enough  to  damp  out  the  oscillations  caused 
by  the  dynamic  forces,  the  flow  is  laminar  and  the  Reynolds  number  is  small. 
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Conversely,  the  flow  is  turbulent  when  the  dynamic  forces  overcome  the  viscous  forces. 
Also,  the  Reynolds  number  is  large  for  turbulent  flow.  [Ref.  22:p.  212] 

The  velocity  profiles  within  the  laminar  and  turbulent  boundary  layers  show  some 
significant  differences.  The  magnitude  of  velocity  in  the  turbulent  boundary  layer  is 
notably  greater,  especially  near  the  reentry  body's  surface.  This  implies  that  the 
turbulent  boundary  layer  yields  much  more  energy  near  the  surface  than  does  a 
corresponding  laminar  flow.  The  velocity  profile  determines  the  skin  friction  on  the 
surface  and  it  can  be  expected  that  greater  heat  transfer  will  occur  under  conditions  of 
a  turbulent  boundary  layer  [Ref.  22 :p.  213].  Figure  39  shows  the  boundary  layer  types. 
Figure  40  shows  boundary  layer  velocity-distance  profiles.  Figure  41  shows  the  altitude, 
air  speed  and  dynamic  pressure  [Ref.  22:pp.  213-214]. 

A  Reynolds  stress  turbulent  boundary  layer  model  which  specifically  accounts  for 
surface  roughness  effects  is  described  by  reference  [84].  In  this  study,  surface  roughness 
is  represented  by  distributed  "sources"  and  "sinks"  in  the  various  governing  equations. 
The  most  significant  term  is  a  sink  term  in  the  mean  momentum  equation,  which 
represents  "form  drag"  on  the  roughness  elements.  [Ref.  84:p.  2] 

A  fundamental  assumption  of  this  model  is  that  the  flow  around  the  individual 
roughness  elements  (only  distributed  roughness  is  considered)  is  attached  to  the  elements. 
[Ref.  84:p.  3]  The  roughness  elements  provide  a  distributed  sink,  due  to  drag,  for  the 
momentum  equation  and  distributed  sources  for  mean  turbulent  kinetic  energy  and 
dissipation.  This  model  also  assumes  that  the  roughness  elements  occupy  no  volume; 
therefore,  this  assumption  becomes  more  severe  as  the  roughness  density  increases.   In 
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[Ref.  22] 
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order  to  compensate  for  this,  the  model  has  been  extended  to  account  for  the  blockage 
effect  of  the  roughness  elements.  [Ref.  84:p.  4] 

The  mean  momentum  equation  is  given  as  [Ref.  84:pp.  4-5] 


ax  ay  ax    ay 


dU 
dy 


by  '    2  Dp 


\-lPT 


(151) 


4/2 
The  variables  U,  V,  u/  and  v/  are  the  reduced  variables,  under  the  boundary  layer 

approximation,  as  given  in  reference  [85].  These  variables  and  this  equation  will  not  be 

derived  here  for  the  sake  of  brevity. 

The  major  advantages  of  this  model  are:  [Ref.  84:p.  5] 

1.  Solutions  are  obtained  for  both  velocity  and  thermal  variables. 

2.  Heat  transfer  is  obtained  directly,  without  invoking  Reynolds  analogy. 

3.  Finite  difference  solutions  are  obtained  using  the  boundary  conditions  that 
fluctuating  quantities  are  zero  at  the  solid  wall  and  in  the  free  stream. 

Reference  [84]  concludes  with  a  comparison  of  a  smooth  wall  turbulent  boundary 
layer  model  and  the  developed  rough  wall  model.  The  rough  wall  model  is  determined 
to  show  that  roughness  spacing  is  more  critical  than  roughness  height,  under  the 
conditions  tested;  however,  the  limited  skin  friction  data  obtained  in  the  study  cannot  be 
interpreted  unambiguously.  [Ref.  84 :p.  38] 
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E.      REENTRY  HEATING  EQUATIONS 

In  Chapter  II,  the  fundamentals  of  reentry  heating  were  discussed  as  presented  by 
the  authors  of  references  [20],  [40]  and  [43].  Aerodynamic  heating,  as  it  applies  to  the 
reentry  of  space  vehicles,  takes  its  roots  in  the  work  of  Allen  and  Eggers,  reference  [65]. 
Reentry  heating  becomes  an  important  consideration  of  the  overall  reentry  process  for 
the  following  reasons:  [Ref.  42:p.  139] 


1 .  Structural  performance  of  the  reentry  vehicle  is  dependent  upon  the  dynamic 
pressures  encountered  during  the  reentry,  which  is  a  function  of  the  reentry 
trajectory. 

2.  Structural  strength  of  the  vehicle  is  a  function  of  the  stresses  induced  by 
temperature  gradients  within  the  component  materials. 

3.  Temperature  gradients  are  proportional  to  the  time  rate  of  heat  input  and 
maximum  time  rate  of  heat  input. 


Therefore,  the  three  critical  parameters  of  the  reentry  trajectory  are  the  total  heat  input 
along  the  trajectory,  the  maximum  rate  of  aerodynamic  heating  and  the  maximum 
dynamic  pressure.  [Ref.  42:p.  139] 

The  mechanism  of  heat  flow  into  the  reentry  vehicle  during  atmospheric  entry  was 
first  described  in  reference  [65],  Since  then,  numerous  combinations  of  reentry  speed 
regimes  and  aerodynamic  shapes  have  led  to  the  publication  of  numerous  technical 
reports.  However,  the  basic  aspects  of  aerodynamic  heating  during  reentry  are  still  the 
same.  The  numerical  factors  for  different  heat  transfer  formulas  and  their  ranges  of 
validity  in  terms  of  the  regime  of  speed  are  the  only  variation  among  the  numerous 
authors.  [Ref.  42:p.  139] 
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The  basic  equations  of  reentry  heating,  equations  (18),  (19)  and  (20),  are  developed 
with  their  simplifying  assumptions  as  presented  in  Chapter  II,  pp.  46-48.  These 
assumptions  yield  the  following  limitations:  [Ref.  42:pp.  141-142] 


1.  Neglecting  radiative  heat  transfer  from  the  vehicle  or  to  the  vehicle  from  the 
high-temperature  air  between  the  shock  wave  and  the  vehicle  surface,  is  based 
on  the  fact  that  the  maximum  allowable  surface  temperature  is  about  the  same 
for  a  variety  of  reentry  vehicle  shapes.  Thus,  outward  radiation  from  the 
surface  will  be  about  the  same.  Neglecting  the  radiative  heat  transfer  from  the 
disturbed  air  is  a  qualitative  simplification  and  therefore  negates  the  application 
of  the  equations  to  very  blunt  and  heavy  shapes  at  reentry  speeds  of  3  km/sec 
or  greater. 

2.  Neglecting  the  real-gas  effects  in  the  flow,  most  importantly  dissociation,  on 
convective  heat  transfer  is  a  good  approximation  for  reentry  speeds  of  up  to  3 
km/sec.  Nevertheless,  this  assumption  is  conservative  and  results  in  higher 
calculated  heating  rates  than  actual  rates. 

3.  Neglecting  the  shock- wave  boundary-layer  interactions  implies  that  the  laminar 
skin-friction  coefficient,  on  a  flat  plate  at  zero  incidence,  is  being  held  constant. 
This  assumption  is  not  valid  at  reentry  speeds  over  6  km/sec. 

4.  Assuming  a  Reynolds  analogy  and  holding  the  Prandtl  number  constant  also 
restricts  the  validity  of  the  equations  to  reentry  speeds  of  less  than  3  km/sec. 


For  the  case  of  low  earth  orbits  and  naturally  decaying  satellites,  these  assumptions 
create  severe  restrictions.  The  circular  orbital  velocity  may  be  calculated  from  [Ref. 
17:p.  38] 


V   = 

c 


M  (152) 

r 


145 


where 

r        =  radius  from  center  of  Earth 

Vc      =  circular  velocity 

Using  this  equation,  it  is  possible  to  approximate  the  circular  orbital  radius  of  a 
satellite  traveling  at  an  altitude  of  near  120  km  or  radius  of  6499  km.  This  is  the  altitude 
of  concern  for  the  focus  of  this  thesis  and  it  can  be  shown  that  the  orbital  velocity  of 
interest  is  approximately  7.8  km/sec.  The  orbital  radius  of  a  satellite  travelling  at  speeds 
of  3  km/sec  and  6  km/sec  are  44,289  km  and  11,072  km  respectively.  This  equates  to 
circular  orbital  altitudes  of  37,910  km  and  4,694  km  . 

The  assumptions  of  reference  [65]  which  are  maintained  in  references  [14],  [40] 
and  [43]  should  be  removed  for  an  accurate  quantitative  analysis  of  the  aerodynamic 
heating  during  the  reentry  of  a  specific  vehicle.  [Ref.  42 :p.  142] 

F.   STRUCTURAL  BREAKUP  OF  A  REENTRY  BODY 

The  previous  sections  of  this  chapter  have  served  to  develop  the  "  state-of-the-art, " 
as  determined  through  the  literature  survey,  of  reentry  formulations  and  solutions.  The 
culmination  of  the  uncontrolled  reentry  process  is  usually  the  structural  breakup  of  the 
reentry  body.  The  breakup  of  a  reentry  body  is  stated  to  be  a  function  of  surface 
temperature,  in  that,  structural  failure  (breakup)  is  assumed  or  expected  to  occur  when 
the  outer  structure  reaches  its  melting  temperature.  [Refs.  63,  86]  It  has  been  shown 
previously  that  the  maximum  temperature  a  reentry  body  will  experience  is  determined 
by  the  maximum  heating  rate,  which  is  a  function  of  ballistic  coefficient.  [Ref.  86:p.  C- 
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11]  Although  the  peak  heating  rate  increases  as  the  ballistic  coefficient  increases, 
average  heating  rate  (a  measure  of  overall  survivability)  and  maximum  local  stagnation 
region  heating  (a  measure  of  local  "hot  spots"),  must  be  considered  in  determining  where 
breakup  occurs.  [Ref.  40:pp.  181-182] 

In  the  1970's,  the  Air  Force  conducted  reentry  experiments  that  used  optical  and 
radar  techniques  to  observe  actual  breakup  events.  The  objective  of  these  tests  was  to 
determine  the  survivability  of  reentry  body  debris.  Specific  findings  from  these  tests 
indicate:  [Ref.  86:pp.  C-4--C-10] 


1.  Classical  convective  heat  transfer  analysis  underestimates  the  reentry  body 
survivability.  Specifically:  actual  breakup  is  at  an  altitude  of  at  least  10  nm 
lower  than  predicted;  actual  surface  temperature  is  at  an  altitude  10  nm  lower 
than  predicted;  and  the  effective  heating  rate  input  is  a  factor  of  four  lower  than 
predicted. 

2.  Consistent  catastrophic  failure  of  magnesium/aluminum  structures  is  at  an 
altitude  of  42  nm.  Figure  42  [Ref.  86]  shows  convergence  of  three  ballistic 
coefficient  lines  close  to  the  magnesium/aluminum  melt  zone  at  this  approximate 
altitude. 

3.  Phenomenon  of  breakup  process  is  independent  of:  body  attitude  and  rates;  body 
diameter;  body  shape;  and  entry  flight  path  angle  (0.3°  >  1.5°). 

4.  Ballistic  coefficient  and  material  of  a  body  determines  survivability. 

5.  A  body  with  low  ballistic  coefficient  will  survive  reentry,  as  will  a  body  with 
a  higher  melting  temperature  survive  reentry  with  a  higher  ballistic  coefficient. 

6.  Reentry  body  structural  integrity  is  maintained  until  melting  temperature  is 
achieved. 

7.  Surface  structure  temperature  is  determined  by  radiation  equilibrium  that  is 
based  on  a  shallow  path  angle  and  a  low  thermal  capacity  of  the  outer  structure. 
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Figure  42:  Structural  Breakup  -  Heating  Rate  vs  Altitude 
[Ref.  86] 


148 


Reference  [85]  concludes  that  classical  convective  heat  transfer  analysis  in  a  free- 
molecular  flow  regime  is  not  indicative  of  the  transition/continuum  flow  heating  which 
is  responsible  for  structural  breakup.  This  finding  is  in  agreement  with  reference  [42] 
which,  as  shown  previously,  discounts  the  classical  reentry  heating  equations  based  on 
the  limitations  imposed  by  the  simplifying  assumptions.  [Ref.  42:pp.  141-142]  Reference 
[64]  also  cites  the  1970's  Air  Force  studies  and  indicates  the  deficiencies  of  the  classical 
convective  heat  transfer  analysis. 

The  consistent  underestimation  of  survivability  was  further  investigated  in  reference 
[64].  By  using  a  calibrated  heating  model,  the  observed  breakup  altitude  of 
approximately  42  nm  (79.5  km)  was  successfully  predicted  for  satellites  with  aluminum 
structures.  Figure  43  [Ref.  64:p.  48]  shows  the  predicted  breakup  altitudes  for  two 
satellite  reentry  simulations.  [Ref.  64:pp.  7-13] 

The  heating  model  was  interfaced  with  a  trajectory  simulation  model  in  order  to 
estimate  the  altitudes  where  possible  breakup  events  could  occur.  Table  XII  shows  the 
breakup  analysis  results  of  the  critical  elements  with  their  associated  material 
composition,  heating  rates  and  breakup  altitudes   [Ref.  64:p.59]. 

Since  the  breakup  analysis  was  assumed  with  strictly  aerodynamic  heating  for 
structures  with  low  melting  points,  the  authors  recommend  further  investigation  for  the 
combined  effect  of  both  aerodynamic  heating  and  aerodynamic  loading  on  breakup. 
Specifically,  for  the  case  where  higher  melting  point  materials  will  survive  to  altitudes 
where  deceleration  loads  are  significant.  [Ref.  64:p.  21] 
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Table  XII: 

[Ref.  64] 


RESULTS  OF  BREAKUP  ANALYSIS 


Configuration 

Critical 

Candidatt 

W* 

\a* 

Breakup 

Eleacnt 

Cooi true t ion 

Altitude 

KactrUl(i) 

(8TU/(c:-») 

(BTu7fc2-a) 

(run.) 

Sp*c*cr*ft 

Boo* 

(DFitargUai 

Epory 

1 

12  -  18 

55 

(2)Hatal  Katria 

8 

12  -  18 

43 

lUaetor/Htat 

Heat  radiacor 

(l)Mo 

688 

69 

- 

Ladiator 

(2)Ti 

142 

69 

- 

(3)Bt 

58 

69 

28 

(4)C 

~ 

69 

- 

Reactor/Radiation 

Radiation 

Ft 

$4 

80  -  104 

31-32 

Shit  Id 

Shit  Id 

Prcnun  Vtntl 

Will 

(l)Mb 

282 

165  -  255 

- 

(2)JUtar 

~ 

• 

~ 

Futl  Pin  #1 

Wall 

Mb 

282 

115  -  1410 

a 

Fuel  Pin  #2 

Wall 

Mb 

282 

247  -  1185 

b 

a        Breakup  anticipated  for  tumbling  mode,  some 
pins  may  survive. 

b        Breakup  anticipated  in  tumbling  mode. 
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Another  series  of  investigations  of  the  breakup  process  are  references  [10]  and  [63]. 
In  these  studies,  the  reentry/breakup  of  Skylab  was  reconstructed  by  piecing  together  all 
available  data  -  after  the  fact.  Both  of  these  reports  show  conclusively,  through  telemetry 
analysis,  that  the  survivability  of  Skylab  was  underestimated.  Reference  [10]  cites  initial 
breakup  predictions  of  120  km  and  shows  that  breakup  did  not  occur  until  at  least  100 
km.  References  [63]  and  [10]  disagree  on  exactly  where  the  OWS  SAS  (Orbital 
Workshop  Solar  Array  System)  separated  from  the  main  body,  based  on  telemetry 
received  at  Ascension  Island.  Reference  [10]  states  that  the  array  was  intact  over 
Ascension  Island;  however,  it  was  not  generating  its  predicted  output.  Therefore,  it  was 
concluded  that  the  array  was  either  bent  back  or  physically  damaged. 

Reference  [63]  concludes  that  the  OWS  SAS  was  completely  separated  from  the 
main  body  prior  to  telemetry  acquisition  at  Ascension  Island  [Ref.  63:p.  344].  Finally, 
reference  [63]  postulates  a  probable  breakup  scenario  as  follows:  [Ref.  63:p.  344] 

1.  OWS  SAS  array  (aerodynamically)  off  at  62  nm  /  117  km. 

2.  ATM  separates  from  remaining  OWS  at  54  nm  /  102  km. 

3.  ATM  SAS  arrays  separate  between  50  and  54  nm  /  94.5  and  102  km. 

4.  ATM  and  OWS  breakup  at  42  nm  /  77.8  km. 

Reference  [10]  concludes  that  breakup  did  not  start  until  some  altitude  below  100  km. 
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IV.     DETERMINISTIC  REENTRY/IMPACT  PREDICTION  METHODS 

A.      CURRENT  REENTRY/IMPACT  PREDICTION  METHOD 

The  following  section  of  this  chapter  describes  the  various  techniques  or  methods 
used  by  different  countries  or  organizations  dealing  with  reentry  and  impact  prediction 
of  naturally  decaying  objects.  It  is  therefore  useful  to  the  reader  to  understand  the 
current  U.S.  method  used  in  reentry /impact  prediction,  since  it  is  the  standard  by  which 
all  other  methods  (as  determined  by  this  literature  survey)  are  compared. 

As  stated  in  previous  chapters,  the  current  method  for  predicting  reentry  time  and 
impact  location  is  that  of  the  Space  Surveillance  Center  (SSC)  located  at  Cheyenne 
Mountain  Air  Force  Base,  or  commonly  referred  to  in  the  literature  as  NORAD.  [Ref. 
87]  NORAD  produces  "element  sets"  which  are  mean  values  of  the  orbital  elements  that 
have  been  obtained  by  removing  the  periodic  orbital  variations  in  a  particular  manner. 
In  order  to  use  these  element  sets,  and  obtain  reasonable  predictions,  these  periodic 
variations  must  be  "reconstructed"  by  the  prediction  model  in  precisely  the  same  manner 
as  they  were  removed  by  NORAD.  Therefore,  an  input  of  NORAD  element  sets  into 
another  model  (even  though  it  may  be  more  accurate,  or  even  into  a  numerical  integrator) 
will  result  in  degraded  predictions  unless,  as  previously  stated,  the  new  model  can 
"reconstruct"  the  periodic  variations.  [Ref.  87:p.  1] 

NORAD  element  sets  are  generated  with  a  general  perturbations  (GP)  model  called 
SGP4.    SGP4  was  developed  by  Ken  Cranford  in  1970.    This  model  was  a  result  of 
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simplification  of  the  more  extensive  analytical  theory  of  Lane  and  Cranford  (1969)  which 
uses  the  solution  of  Brouwer  (1959)  for  its  gravitational  model  and  a  power  density 
function  for  its  atmosphere  model.  [Ref.  87:p.  3]  It  should  be  emphasized  that  this 
atmosphere  "model"  is  a  static  representation  of  density  as  opposed  to  the  dynamical 
models  discussed  in  Chapter  II,  such  as  those  of  Jacchia  and  others. 

The  gravitational  model  includes  J2  and  J3  harmonics;  however,  J4  and  J5  were 
dropped  in  order  to  avoid  singularities  occurring  at  critical  inclination  [Ref.  88 :p.  2]. 
Rates  of  change  of  mean  motion  and  eccentricity  are  derived  from  the  density  function. 
The  product  of  ballistic  coefficient  and  a  reference  density,  denoted  B*,  is  treated  as  a 
solved  for  parameter.  Coupling  between  J2  and  drag  is  included  in  the  argument  of 
perigee,  right  ascension  of  the  node  and  mean  anomaly.  The  mean  motion  of  SGP4  is 
a  pure  Brouwer,  or  two-body,  mean  motion.  [Ref.  88:p.  3] 

Reentry  and  impact  predictions,  in  the  U.S.,  are  made  using  a  special  perturbations 
(SP)  propagator  with  conversions  between  GP  and  SP  theories  handled  as  outlined  in 
references  [87,  89-90].  The  GP  theory  is  fast,  analytical  and  of  low-accuracy,  when 
compared  to  the  SP  theory.  SP  theory  uses  a  Gauss- Jackson,  eighth-order,  numerical 
integrator,  incorporates  a  6- 12th  order  geopotential  model  and  applies  a  dynamical 
atmospheric  density  model  (Jacchia-65).  The  "conversion"  between  GP  and  SP  theories 
is  the  process  of  performing  an  SP  differential  correction  of  the  initial  state  vector  as 
derived  from  the  GP  theory  (NORAD  two-line  element  set).  These  are  the  initial 
standards  for  GP  and  SP  theory  compatibility  which  must  be  considered  in  the  following 
discussion  of  different  reentry /impact  prediction  methods.  [Ref.  91:p.  8] 
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The  first  TIP  (Tracking  and  Impact  Prediction)  run  is  performed  approximately  15- 
20  days  prior  to  the  estimated  reentry  date,  which  is  initially  predicted  by  the  GP  model. 
This  is  also  when  tasking  of  the  observation  sites  is  initially  increased  in  order  to  support 
high  accuracy  SP  processing.  [Ref.  92 :p.  1] 

Orbit  determination  is  accomplished  through  a  first-order,  linear,  weighted,  least 
square,  curve  fitting  process,  commonly  called  differential  correction.  Sliding  fits  are 
used  to  process  both  new  and  old  metric  data  until  the  satellite  is  "no  longer  in  orbit." 
The  force  models  used  are  the  Jacchia  dynamic  atmosphere  (1965)  and  the  World 
Geodetic  System  Earth  gravity  model  (1972).  The  Earth  gravity  model  is  truncated  to 
the  sixth  order  for  satellite  decay  predictions.  [Ref.  92:p.  2] 

Direct  step-by-step  numerical  integration  of  the  total  acceleration  acting  on  the 
decaying  satellite  is  accomplished  in  the  manner  of  Cowell's  method.  Gauss- Jackson 
eighth  order  predictor-corrector  formulas,  in  ordinate  form,  are  used  to  integrate  the 
equations  of  motion.  Because  of  computer  run-time  constraints,  the  partial  derivatives 
necessary  for  differential  correction  are  computed  analytically  except  for  the  secular 
variations  due  to  atmospheric  drag  of  the  orbit  semi-major  axis  and  eccentricity.  These 
parameters  are  integrated  numerically  using  a  low-order  (trapezoidal  rule)  integrator. 
[Ref.  92:p.  2] 
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Currently,  the  differential  correction  solution  state  consists  of  the  equinoctial 
elements  and  satellite  ballistic  coefficient  model  parameter,  B0.  After  each  solution,  the 
new  state  is  used  to  predict  a  decay  time  (when  altitude  equals  10  km).  [Ref.  92:p.  2] 

B.      ALTERNATE  REENTRY/IMPACT  PREDICTION  METHODS 

1.      Reentry  Prediction  Methods  At  ESOC 

One  method  used  at  the  European  Space  Operations  Center  (ESOC)  for  the 
prediction  of  reentry  and  impact  of  decaying  satellites  is  an  improved  and  computerized 
version  of  the  King-Hele  technique,  reference  [55].  [Ref.  93]  Another  method  used  at 
ESOC  is  based  in  a  computer  program  called  FOCUS  [Ref.  94]. 

The  principal  characteristic  of  FOCUS  is  its  ability  to  overcome  the 
deficiencies  of  the  semi-analytical  orbit  prediction  techniques  at  low  altitudes  [Ref.  94:p. 
26].  Recall  that  in  near-Earth  orbits,  200-700  km,  Earth  oblateness  (J2)  is  regarded  as 
the  only  first-order  perturbation.  Higher  zonal  harmonics  and  air  drag  are  regarded  as 
second-order  contributions.  All  other  effects  are  considered  as  less  than  second-order 
perturbations  (less  than  J22  in  magnitude).  [Ref.  94:p.  25]  Eventually,  the  satellite  passes 
through  an  altitude  regime  where  the  air  drag  force  (caused  by  increasing  atmospheric 
density)  reaches  a  magnitude  of  the  same  order  as  the  J2  Earth  oblateness  effect.  At  this 
altitude,  «  150  km,  depending  on  the  vehicle's  ballistic  coefficient,  the  accuracy  of 
analytically  derived  drag  perturbation  results  strongly  deteriorates.  [Ref.  94:p.  26] 

FOCUS  stops  the  state  propagation  after  passing  through  a  user-defined 
altitude  shell  (h  =  170  km  specifically  for  the  case  of  Salyut-7/Cosmos-1686,  which  is 
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the  focus  of  reference  [94])  and  forwards  a  calculated  osculating  Keplerian  state  vector 
at  epoch,  together  with  all  relevant  perturbation  parameters.  This  new  set  of  data  is  then 
numerically  integrated  and  the  reentry  trajectory  is  propagated  until  shortly  before  impact 
with  the  Earth.  [Ref.  94:p.  26] 

Several  key  features  of  the  FOCUS  program  are:  [Ref.  94:p.  26] 

1.  Perturbation  equations: 

(a)  Cowell's  formulation  of  the  perturbed  Newtonian  equations  written  in  terms 
of  six  first-order,  differential  equations  for  each  component  of  the  Cartesian 
state  vector 

(b)  Reference  frame  is  the  mean  equatorial  system  of  date 

2.  Perturbation  models: 

(a)  Geopotential  model  GEM  10B  (J2-J7  used) 

(b)  Atmospheric  density  models:  MSIS-86  for  altitudes  >  120  km,  U.S. 
Standard  Atmosphere  (USSA-76)  for  altitudes  <  90  km,  and  a  bridging 
function  for  altitudes  between  90  and  120  km. 

(c)  Variable  drag  coefficient,  CD  =  f(Ma,  Re,  Kn) 

(d)  Co-rotating  Earth  atmosphere 

(f)  Luni-solar  third  body  attraction  (point  mass) 

(g)  Solar  radiation  pressure 

3.  Integrator: 

(a)  Runge-Kutta/Shanks  7/8  single  step  method  for  generation  of  a  starting  arc 

(b)  Adams-Bashforth/Adams-Moulton(AB/AM)forth-orderpredictor-corrector, 
multi-step  method  for  propagation  of  the  Cartesian  state  vector. 
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(c)   Non-regularized  time  (t)  used  as  an  integration  variable,  with  constant  step 
sizes  of  At  =  30  sec. 


A  significant  limiting  feature  of  this  program  is  that  the  reentry  trajectory  is 
terminated  at  30  km  altitude  because  the  governing  laws  of  perturbed  Keplerian  motion 
become  invalid  below  this  threshold  altitude.  This  criterion  is  marked  by  a  decrease  in 
the  orbital  energy  to  a  level,  where  the  aerodynamic  forces  are  essentially  in  balance  with 
the  zero-th  order  central  gravitational  attraction  term.  [Ref.  94:p.  26] 

The  reentry  vehicle  is,  however,  considered  to  be  in  nearly  vertical  fall  from 
an  altitude  of  30  km  and  below.  Thus,  there  is  only  a  minor  dispersion  of  the  impact 
point  during  the  final  seconds  of  flight.  This  rationale  leads  to  the  conclusion  that  it  is 
not  necessary  to  perform  another  transition  from  strongly  perturbed  Keplerian  motion  to 
an  aerodynamic  flight  phase  for  the  integration  to  Earth  impact.  Thus  the  Center  Of 
Impact  Window  (COIW)  is  defined  as  the  location  at  which  the  vehicle  passes  through 
the  30  km  altitude.  [Ref.  94 :p.  26] 

The  aerodynamic  transition  regime  is  defined  as  a  computed,  weighted  mean 
of  Pmsis  and  pUSSA  where 

Pmsis  =f(h^A,UTtxttdiFl(lvFl0VAp)  (153) 

Pussa  =/W  <154> 

which  results  in 
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P  =  wmsis  Pmsis  +  0  -w«ffls)  Puss*  (155) 

where  the  weighting  factor 

*W(*>  e  [0,1]  (156) 

is  defined  for  the  altitude  region  h120  ^  h  >  h^  as 


1  (        ^120  "^ 

w*sff(*)  "  tCos|ti 


2      I     /z120  -  h^ 


+ 1  (157) 

2 


where 

h  =  geodetic  altitude 

$  =  geodetic  latitude 

X  =  geographic  longitude 

t  =  local  solar  time 

UT  =  universal  time 

td  =  day  of  the  year 

When  the  altitude  is  high  enough  to  maintain 


Ah.  =  —  >  1  (158) 

dslc 


where 

d,/c     =  characteristic  length  of  vehicle 
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then  CD  is  equal  to  a  constant  given  by  Co101,  provided  Kn„  >  1.  As  the  altitude 
decreases,  a  transition  to  hypersonic  continuum  flow  is  entered,  where  Kn^  =  f(h)  and 
[Ref.  94:p.  27] 

CD  =  C?  +  C?loge(C?KnJ  (159) 

provided  that  0.02  <  Kn^  <  1.  During  the  hypersonic  continuum  phase  another 
constant  CD  level  is  attained  (about  50%  less  than  the  free  molecular  flow  value),  which 
is  given  by 

CD  *  constant  =  C0Ma  <16°) 

provided 

Knw  <  0.02 

Maw  >  5 

The  next  phase  of  supersonic   and   transonic   continuum   flow  can  be 
approximated  by  an  altitude  dependent  Mach  number,  Ma  =  f(h),  where 

CD  =  Cffl  +  C2Wj(Maa>-0.4)**C3Wj  x  exp[-Cf(Ma..-0.4)]  (161) 

provided 

Kn„  >  0.02 
0.8  <  Ma,,  <  5 

Finally,  in  the  subsonic  phase,  the  drag  coefficient  is  dependent  upon  viscous 
interactions  and  therefore  Reynolds  number,  Re  =  f(h),  given  by 
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provided 

Kn,,  <  0.02 
Ma,,  <  0.8 

These  dependencies,  CD  =  f(Kn,  Re,  Ma),  with  the  underlying  altitude 
functions  Kn(h),  Ma(h)  and  Re(h)  are  from  the  U.S.  Standard  Atmosphere  (USSA)  and 
are  incorporated  into  the  numerical  reentry  prediction  software.  These  model  constants 
are  limited,  however,  to  spherical  and  cylindrical  shapes,  with  their  longitudinal  axis 
perpendicular  to  the  airflow.  [Ref.  94 :p.  27] 

Since  Ma  =  V,to  /  V^,,  and  Re  =  Vatm  d/v  are  functions  of  geodetic  latitude 
related  by  V^^n),  kinematic  viscosity,  v(h),  and  aerodynamic  velocity,  V,^,  the  free 
fall  of  the  reentering  vehicle  in  the  lower  atmosphere  may  be  determined  by  iteratively 
solving  the  following  equation  for  the  equilibrium  descent  velocity,  Vltm 

li-*a-25 1© (163) 

~  Ip(A)cd(UJ 

where 

g(h)    =  central  gravitational  acceleration  of  the  Earth 

p(h)    =  local  air  density 

Using  this  technique,  in  the  post-flight  analysis  of  Salyut-7/Cosmos-1686,  the 
following  values  are  determined: 


161 


m/A   =  159.5  kg/m2 
therefore  V,,,,,  at 

h        =30  km,  is  equal  to  330  m/s 

h        =20  km,  is  equal  to  100  m/s 

h        =10  km,  is  equal  to  70  m/s 

Figure  44  shows  the  ESOC  reconstruction  of  the  Salyut-7/Cosmos-1686  final 
descent  altitude  profile.  The  actual  reentry  occurred  at  0345  UCT  on  7  Feb.  1991  over 
(39.3°S,  69.7°W)  [Ref.  94:p.  29].  Table  13  shows  a  comparison  of  reentry  predictions 
from  the  U.S.  and  various  European  and  former  Soviet  sites  [Ref.  94 :p.  32]. 

The  conclusions  of  this  European  Space  Agency  (ESA)/ESOC  report  indicate 
that  there  was  good  correlation  between  ESA,  U.S.  and  former  Soviet  predictions.  There 
was  some  difficulty  in  interpreting  the  Soviet  orbit  determination  results  at  approximately 
12  hours  prior  to  reentry.  Also  at  this  time,  the  U.S.  element  sets  became  "time  late," 
due  to  transmission  delay  times  from  the  U.S.  to  ESOC.  Therefore,  the  ESA/ESOC 
predictions  could  not  be  maintained  in  "real  time"  after  this  point.  [Ref.  94:p.  33] 

2.       The  LIFETIME  Model 

The  Aerospace  Corporation  previously  developed  the  LIFETIME  program 
and  it  was  further  refined  through  the  work  of  reference  [95].  This  program  is  similar 
in  some  respects  to  the  previously  discussed  ESA  program,  FOCUS.  This  program 
offers  a  "fast,  efficient  computer  tool  for  orbital  lifetime  estimation"  [Ref.  95:p.  17]. 
An  advertised  major  benefit  of  the  LIFETIME  program  is  its  usefulness  as  a  "highly 
accurate,  real-time  reentry  prediction  tool."  [Ref.  95 :p.  20]     This  may  also  have 
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Hgure  44:  Salyut-7/Cosmos-1686  Final  Descent  Altitude  Profile 
[Ref.  94] 
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Table  XIII:  SALYUT-7/COSMOS-1686  FINAL  DECAY  PREDICTIONS 
[Ref.  94] 


Prediction 
Source 

Re- Entry  Predictions 
(times  in  UTC) 

Received 
at  ESOC 

Re-Entry  Pred. 
(COIW) 

NASA  GSFC 

01-Feb  14:58 

06-Feb  21:33 

ESOC  MAS 

02-Feb  04:27 

06-Feb  22:30 

CNES  CST 

03-Feb  11:33 

06-Feb  04:00 

ESOC  MAS 

04-Feb  10:30 

07-Feb  01:37 

NASA  CSFC 

04- Feb  16:04 

07-Feb  03:38 

CNUCE  CNR 

04-Feb  19:43 

07-Feb  03:17 

CNES  CST 

05-Feb  04:27 

06-Feb  18:23 

ESOC  MAS 

05-Feb  10:00 

07-F*h  02:23 

NASA  CSFC 

05-F«b  15:05 

07-Feb  03:05 

ESOC  MAS 

06-Feb  11:33 

07-Feb  03:50 

RAE 

06-Feb  11:38 

07-Feb  04:30 

CNUCE  CNR 

06-Feb  12:50 

06-Feb  22:26 

NASA  CSFC 

06-Feb  14:09 

07-Feb  03:28 

CNUCE  CNR 

06-Feb  16:51 

07-Feb  01:24 

NASA  CSFC 

06-Feb  16:58 

07-Feb  04:14 

CNES  CST 

06-Feb  17:00 

07-Feb  05:00 

ESOC  MAS 

06-Feb  19:47 

06-Feb  04  19 

CNUCE  CNR 

06-Feb  22:05 

07-Feb  04:43 

CNES  CST 

06-Feb  22:15 

07-Feb  04:38 

IKI  IAM 

06-Feb  23:37 

07-Feb  04:40 

NASA  CSFC 

06-Feb  23:44 

07-Feb  04:05 

ESOC  MAS 

07-Feb  0100 

07-Feb  04  50 

IKI  IAM 

07-Feb  01:02 

07-Feb  04:38 

IKI  MCC 

07-Feb  01:02 

07-Feb  04:20 

CNUCE  CNR 

07-Feb  01:07 

07-Feb  04:00 

IKI  MoD 

07-Feb  01:57 

07-Feb  04:00 

IKI  MCC 

07-Feb  02:22 

07-Feb  03:57 

IKI  IAM 

07-Feb  04:00 

07-Feb  03:37 

ESOC  MAS 

07-Feb  04  20 

07-Feb  04:29 

IKI  MoD  * 

07-Feb  04:34 

07-Feb  03.47 

ESOC  MAS  ■ 

07-Feb  13  00 

07-Feb  03:45 
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application  as  a  first  time  modeling  tool  as  compared  to  methods  capable  only  of  post 
flight  analysis  of  reentry  events. 

The  LIFETIME  program  uses  either  the  Jacchia-Walker  (64)  or  Jacchia  (71) 
atmospheric  density  models  for  altitudes  above  90  km  and  uses  the  U.S.  Standard 
Atmosphere  (1962)  for  altitudes  below  90  km.  The  latest  version  (4.0),  the  subject  of 
reference  [93],  also  allows  differential  corrections  of  the  ballistic  coefficient  and 
movement  of  solar  panels  to  simulate  sun  tracking.  [Ref.  95:pp.  18-19]  Version  4.0  was 
designed  to  solve  a  major  deficiency  of  version  3.0,  in  that,  there  was  a  built-in 
uncertainty  in  Earth  impact  time  prediction  of  at  least  one  orbital  revolution.  This  was 
inherent  due  to  the  minimum  step  size  of  one  revolution,  which  was  a  function  of  using 
the  averaged  equations  of  motion.  [Ref.  95 :p.  22] 

A  "unique"  feature  of  the  LIFETIME  program  is  the  differential  correction 
of  the  ballistic  coefficient.  By  using  the  least  squares  method,  the  following  equation  for 
differentially  correcting  the  inverse  ballistic  coefficient  (B*)  can  be  formulated  [Ref. 
95:pp.  28-29] 


Afl*  = 


^'   {     '  dB*  ldB\ 
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(  de, 


\2 


[dB') 
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where 


N       =  number  of  observations 


=  a,/ao 
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aj        =  observed  semi-major  axis  at  time  i 

ao       =  semi-major  axis  at  epoch 

e{        =  observed  eccentricity  at  time  i 

It  can  be  shown  that  orbit  decay  in  semi-major  axis  and  eccentricity  is 
directly  proportional  to  the  product  of  inverse  ballistic  coefficient  and  density  (B*p). 
This  differential  correction  process  (which  absorbs  the  atmospheric  density  uncertainties) 
results  in  a  "converged  ballistic  coefficient"  through  multiple  iterations  of  equation  175. 
This  is  how  LIFETIME  becomes  more  accurate  than  the  atmosphere  model  it  uses.  [Ref. 
95:p.  29] 

LIFETIME  integrates  the  averaged  equations  of  motion  using  a  step  size  of 
multiples  of  the  orbital  period,  until  the  satellite  reaches  a  specified  or  default  altitude. 
Based  on  [86],  the  default  altitude  chosen  is  42  nm  (78  km)  and  this  is  where  "breakup" 
is  said  to  occur.  LIFETIME  then  backs  up  to  the  last  propagation  step,  resets  the  orbital 
elements  and  time  to  the  previous  step's  values.  Next,  the  Runge-Kutta  7(8)  integration 
routine  is  invoked  and  the  satellite's  orbit  is  propagated  through  breakup  to  impact.  (The 
elements  are  converted  from  classical  mean  to  classical  osculating  elements  and  finally 
to  Cartesian  elements  during  the  numerical  integration  process.)  [Ref.  95:pp.  32-60] 
Numerical  integration  is  the  preferred  method  for  dealing  with  the  final  stages  of  reentry 
to  impact  since  it  is  the  "most  efficient  and  accurate  way  to  handle  regions  where  a  high 
rate  of  change  in  the  equations  of  motion"  are  prevalent  [Ref.  95 :p.  55]. 
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Several  significant  user  input  altitude  values  are:  [Ref.  95 :p.  59] 


1.  RKALT  -  lower  limit  at  which  LIFETIME  converts  to  Cartesian  numerical 
integration.   Default  =  0.000  km,  most  effective  if  >  BRKALT. 

2.  BRKALT  -  vehicle  breakup  altitude.   Default  =  77.784  km. 

3.  END  ALT  -  perigee  altitude  lower  limit,  ends  integration  and  propagation. 
Default  =  10.000  km,  this  should  be  set  to  0.000  km  to  model  impact. 


At  the  point  where  the  satellite  reaches  the  BRKALT,  the  lat/long  projection 
of  the  satellite  is  recorded.  This  point  is  defined  as  the  "debris  heel  point,"  or  the  first 
point  along  the  groundtrack  where  debris  could  potentially  impact.  Propagation  continues 
until  ENDALT,  where  the  lat/long  is  identified  as  the  center  of  mass  impact  point. 
LIFETIME  then  resets  the  Cartesian  elements  to  those  at  the  BRKALT  point,  however, 
B  is  changed  to  60  lb/ ft2  and  the  propagation  to  ENDALT  is  continued.  This  allows  the 
computation  of  a  "debris  toe  point,"  or  the  point  of  farthest  travel  of  any  breakup  debris. 
The  value  of  60  lb/ft2  for  the  ballistic  coefficient  is  based  on  reference  [86].  [Ref.  95:pp. 
59-60]  Figure  45  shows  the  debris  dispersion  footprint  as  described  above  [Ref.  95  :p. 
60].  Figure  46  shows  the  LIFETIME  groundtrack  and  impact  area  plot  as  well  as  a 
sample  altitude  decay  history  [Ref.  95 :p.  63].  (An  item  of  interest  not  addressed  by  the 
author  of  reference  [95]  is  the  fact  that  the  center  of  mass,  in  this  sample  plot,  impacts 
after  the  "final  debris  impact.") 

The  results  of  impact  prediction  accuracy  comparisons  between  LIFETIME 
3.0  and  4.0  do  not  definitively  show  4.0  as  a  significant  improvement  over  3.0. 
However,  some  notable  improvement  in  accuracy  for  version  4.0  can  be  shown  and 
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Figure  45:  LIFETIME  Debris  Dispersion  Footprint 
[Ref.  95] 
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Figure  46:  LIFETIME  Groundtrack  And  Altitude  Decay  History 
[Ref.  95] 
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greater  confidence  is  held  in  version  4.0's  period  calculations  and  final  impact 
predictions.  [Ref.  95:pp.  65-70] 

A  sensitivity  study  and  analysis  was  performed  in  order  to  determine  the: 

1.  Effects  of  NORAD  data  span  on  program  accuracy 

2.  Effects  of  prediction  span  on  program  accuracy 

3.  Effects  of  solar  conditions  during  differential  corrections 

4.  Effects  of  solar  flux  on  overall  program  sensitivity 

The  results  of  the  sensitivity  study  and  analysis  are  as  follows:  [Ref.  95:pp.  89-94] 

1.  NORAD  data  span: 

(a)  Period  calculation  accuracy  -  no  strong  trend  for  sensitivity  noted 

(b)  Impact  prediction  accuracy  -  general  trend  for  less  averaged  impact  time 
error  for  greater  time  spans 

2.  Prediction  span  -  impact  time  errors  show  a  strong  trend  for  increased  accuracy 
with  shorter  prediction  spans. 

3.  Solar  flux  -  the  results  indicated  a  complex  relationship  between  the  estimated 
solar  flux  for  the  prediction  period,  the  actual  values  of  solar  flux,  the 
prediction  span,  and  the  resulting  impact  errors. 

A  general  observation  and  conclusion  is  that  fairly  steady  solar  activity  at 
modest  levels  provide  for  the  most  accurate  impact  prediction,  especially  when  using  the 
solar  flux  value  of  the  last  day  on  NORAD  data  as  a  constant  during  the  prediction  span. 
During  periods  of  highly  variable  solar  flux,  the  last  day's  value  of  NORAD  data  may 
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not  give  the  least  error,  however,  this  conservative  assumption  will  not  result  in  the 

maximum  error  either.  [Ref.  95: p.  92] 

An     appropriate     concluding     remark     of     reference     [95]     is     that 

It  seems  there  are  still  dynamic  aspects  of  the  final  orbit  decay  and  reentry  process 
that  warrant  further  study.  [Ref.  95 :p.  93] 

3.       Modeling  Ballistic  Coefficient 

Analysis  of  historical  data  shows  that  in  the  final  days  prior  to  decay,  the 
ballistic  coefficient,  B,  (CDA/m)  exhibits  a  pronounced  variation  with  time.  This 
variation  is  due  to  one  or  more  of  several  causes:  [Ref.  92 :p.  2] 

1.  Tumbling  or  weather  vaning  of  the  satellite 

2.  Loss  of  mass 

3.  Variation  of  the  drag  coefficient 

Presently,  none  of  these  effects  are  modeled  explicitly  since  this  requires  precise 
knowledge  of  the  multiple  variables  described  in  Chapters  II  and  III,  such  as  vehicle 
mass,  attitude  and  shape.  Additionally,  any  unmodeled  forces  acting  in  the  in-track 
direction,  such  as  those  inaccuracies  in  the  atmospheric  density  model,  can  manifest 
themselves  in  the  solution  ballistic  coefficient.  [Ref.  92 :p.  2] 

Despite  the  lack  of  specific  knowledge  of  unknown  forces  acting  on  a 
decaying  satellite,  an  improvement  in  reentry  prediction  accuracy  can  be  made  by 
introducing  a  new  model  parameter.  This  parameter  takes  into  account  the  fact  that  for 
many  decaying  satellites,  the  variation  in  B  appears  to  be  linear  during  the  last  several 
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days  prior  to  reentry.   This  new  parameter  is  B0-dot,  the  time  derivative  of  the  ballistic 
coefficient,  which  is  parameterized  as 

B  =  BQ  *  BQ(f-t^  (165) 

where  subscript  0  indicates  the  quantities  at  the  solution  state  epoch,  to-  [Ref.  92:p.  3] 
The  introduction  of  B0-dot  requires  that  new  semi-analytic  partial  derivatives 
for  B0  and  B0-dot  be  calculated.  These  must  be  determined  over  a  finite  time  interval 
using  a  differential  correction  fitting  process.  Several  such  differential  corrections  must 
be  made  over  the  course  of  the  final  days  of  orbit  decay  and  reentry,  with  a  new  B0  and 
B0-dot  being  determined  for  each  differential  correction.  It  is  necessary  to  start  the 
fitting  process  with  close  approximations  of  B0  and  B0-dot  for  the  first  differential 
correction.  These  initial  "guesses"  come  from  analysis  of  historical  data.  [Ref.  92:p.  3] 
Reference  [92]  studied  264  reentry  events  over  the  time  span  from  Jan.  1977 
to  Mar.  1979.  All  of  these  objects  were  processed  as  TIP  reentries  and  historical  records 
from  the  SSC  were  used  in  this  investigation.  In  order  to  optimize  the  data,  the  records 
of  156  reentries  were  chosen,  since  these  represented  rocket  bodies  and  payloads  from 
the  former  USSR,  which  presumably  would  have  well-known  physical  attributes.  Table 
14  shows  the  data  distribution  used  throughout  this  investigation.  [Ref.  92:pp.  3-4] 

The  development  of  the  solutions  for  the  partial  derivatives  of  Bq  are  not 
shown  here  for  the  sake  of  brevity,  however,  the  results  are  as  follows  [Ref.  92:pp.  10- 
11] 
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Table  XIV:  DISTRIBUTION  OF  DECAYED  ORBITAL  OBJECTS 
[Ret.  92] 


DISTRIBUTION  OF  264  TIP  LOG  DECAYED  ORBITAL  OBJECTS 

YEAR:  1970 

MONTH:  J       FMAMJ       J       ASOND 

NUMBER:  2 

YEAR;  1977 

MONTH:  J       FMAMJ       J       ASOND 

NUMBER:  8       8       116        158       8       11149       8       7 

YEAR:  1978 

MONTH:  J       FMAMJ       J       ASOND 

NUMBER:  119       9       136       127       8        11138       18 

YEAR:  1979 

MONTH:  J        FMAMJ        J        ASOND 

NUMBER;  6        10     7 

YEAR:  1983 

MONTH:  J       FMAMJ       J       A      S       O      N      D 

NUMBER;  1 


• 

Selected    156 

Soviet    orbitally 

decayed   object 

distribution 

.6.7.8 

•           ROCKET  BODIES 

DIA(m) 

LNG  (m) 

21 

A2        VOSTOK 

(440    kg) 

2.6 

7.5 

89 

A2        VOSTOiC 

(2500   kg) 

2.0 

2.0 

4 

Bl         COSMOS 

(1500   kg) 

1.65 

8.0 

8 

O         SKEAN 

(2000   kg) 

2.0 

6.0 

3 

Dl         PROTON 

(1900   kg) 

3.9 

4.0 

3 

Dl         PROTON 

(4000  kg) 

4.0 

12.0 

3 

F1M      SCARP 

(  700    kg)   • 

2.0     - 

5.0 

1 

F1M     SCARP 

(1500  kg) 

..     2.5 

6.0 

PAYLOADS5 

3 

NAVSAT 

(  875    kg) 

? 

7 

1 

(  325  .  kg) 

1 

(3640   kg) 

1 

(  680  ,.kg) 

8 

MTLSAT 

(  400  '  kg) 

i 

MIL/SCI 

(550*  kg) 

4 

RORSAT* 

(4500   k'g)- 

2 

MIL 

(900    kg) 

1 

(1120  kg) 

1 

-(3I0.Q*.kg). 

9'  The  exact  size  and  shape  of  the  payloads  are  unknown  at  this  point. 
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bJH  =B0a  n['~dt  (166) 

°dBQ        0J0Jt.Be 

BQ^L=BQa  ^['--dt  (167) 

°dB0       0*°KBe 


B0*L--lBo!!*r±adi  (168) 

0  dB0         2    °  a0  K  B 


*f --!*5*-0 


»  J!*  =»  i*  =0  (170) 

°  dBQ        °  3B0 


where 

aj  =  e  cos  (ft  +  «) 

\  =  e  sin  (ft  +  co) 

e  =  time  rate  of  change  of  eccentricity  due  to  drag  perturbations 

a  =  time  rate  of  change  of  semi-major  axis  due  to  drag  perturbations 

L  =  mean  longitude 

X  =  tan  (i/2)  sin  ft 

¥  =  tan  (i/2)  cos  ft 
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Similarly,  the  partial  derivatives  of  Bo-dot  are  given  [Ref.  92:p.  12] 
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Of  the  156  cases  studied  in  the  development  of  the  new  parameter  and  its 
associated  partial  derivatives,  two  case  studies  were  chosen  for  the  application  of  the  new 
method.  These  cases  were  Cosmos-954  and  Cosmos-1402,  both  of  them  being  significant 
due  to  their  potentially  survivable  nuclear  reactors.  [Ref.  92:p.  13] 

The  conclusion  of  this  paper  is  that  the  new  method  is  more  accurate  than  the 
conventional  method  as  evidenced  by  the  comparison  of  the  best  fit  differential  correction 
root  mean  square  (RMS)  of  residuals  for  each  of  the  25  sliding-fit  runs  over  the  last  nine 
days  of  orbit  for  Cosmos-954.  The  same  comparison  is  made  for  Cosmos-1402  with  the 
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only  exception  that  26  vice  25  SP  differential  corrections  were  run.  The  results  of  these 
comparisons  are  shown  in  Figure  47  [Ref.  92: pp.  15,  17]. 

As  stated  previously,  the  historical  data  analysis  of  most  satellites  showed  the 
linear  trend  in  B  only  over  the  last  3-5  days  prior  to  reentry.  However,  one  of  the  two 
specific  cases  chosen  for  application  of  the  new  method,  Cosmos-954,  did  not  exhibit  this 
trend  until  the  final  day  prior  to  reentry.  Table  15  shows  a  comparison  of  decay  time 
error  root  mean  square  over  the  last  five  days  of  orbit  (for  these  two  case  studies)  as  well 
as  a  similar  comparison  after  the  linear  trend  in  B  becomes  apparent  in  both  reentries. 
[Ref.  92:p.  18] 

Another  method  using  a  correction  of  the  ballistic  coefficient  comes  from 

reference  [47].    This  investigation  states  that  during  the  joint  international  effort  to 

predict  the  reentry/impact  of  Salyut-7/Cosmos-1686,  it  was  observed  that  the  "most  exact 

data"  had  consistently  been  presented  by  the  Space  Control  System  (SCS)  of  the  former 

USSR.   This  was  explained  by  the  fact  that: 

...great  attention  has  been  given  to  the  problem  of  determination  and  the  prediction 
of  the  satellite  motion  in  the  atmosphere:  for  the  solution  of  this  problem  they  have 
been  conducting  the  cycle  of  theoretical  and  experimental  studies  during  [the  last] 
20  years.  [Ref.  47:p.  35] 

The  model  used  by  the  SCS  to  describe  satellite  motion  will,  to  a  great 

extent,  determine  the  "completeness"  and  precision  of  the  time  and  area  of  low-orbit 

satellite  reentry  predictions.    Gravitational  field  models  in  spherical  functions  of  36  or 

higher  orders  are  not  necessary,  they  only  serve  to  increase  the  computation  time.   At 
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Figure  47:  Best  Fit  DC  RMS  For  Cosmos-954/1402 
[Ref.  92] 
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Table  XV:  DECAY  TIME  ERROR  RMS 
[Ref.  92] 


DECAY  TIME  ERROR  ROOT  MEAN  SQUARE  (RMS)  OVER  LAST  5  DAYS 

Decay  Error  RMS(hrs) 
Decay  Error  RMS(hrs)    Bo,  BQ  Solution, 
Satellite  Bo  Solution  B^.o.O  for  Prediction    lmprovement(%) 


COSMOS  954  3.51       .  4.11  -17.09 

COSMOS  1402  1.59  1.43  10.06 

OECAY  TIKE  ERROR  ROOT  MEAN  SQUARE  (RMS)  AFTER  B  BECOMES  LINEAR 

Decay  Error  RMS(hrs) 
Decay  Error  RMS(hrs)    Bo,  kQ  Solution, 
Satellite  Bq  Solution  BQ«o.O  for  Prediction    lmprovement(%) 


COSMOS  954  1.22  0.48  60.66 

COSMOS  1402  1.51  0.98  35.10 
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low  altitudes,  two  to  three  days  prior  to  reentry,  a  simple  geopotential  model  is  suitable, 
since  it  is  the  atmospheric  resistance  force  which  is  increasing  by  several  orders  of 
magnitude,  not  gravitational  forces.  [Ref.  47:p.  36] 

Prediction  of  low-orbit  satellite  motion  is  accomplished  by  numerical  and 
semi-analytical  integration  methods.  The  semi-analytical  methods  are  based  on  the 
asymptotic  solutions  given  by  Krilov-Bogoljubov,  Zeipel-Brouwer  averaging  methods. 
[References  given  are  written  in  Russian  and  the  authors  were  unable  to  obtain  English 
translations.]  [Ref.  47:p.  36] 

The  actual  method  applied  to  a  reentry  event  may  be  a  variant  of  some 
numerical  and  semi-analytical  methods.  These  specialized  fast-acting,  semi-analytical 
prediction  algorithms  for  low-altitude  orbits  have  "counting  times  of  order  0.01  seconds" 
for  a  24  hour  prediction  span.  [Ref.  47:p.  36] 

The  reentry  of  the  Salyut-7  complex  was  predicted  and  confirmed  by 
observation  in  the  range  of  75-105  km.  In  this  altitude  range,  it  was  also  observed  and 
predicted  that  breakup  would  occur.  [Ref.  47:p.  39] 

It  was  assumed  that  the  ballistic  coefficient  would  vary  during  reentry,  to  the 
same  degree  as  it  did  prior  to  an  attitude  control  maneuver.  This  assumption  allowed  the 
reentry  problem  to  be  solved  by  varying  the  initial  conditions  within  their  ranges  of 
possible  errors,  with  subsequent  integration  of  motion  equations  up  to  the  point  of 
reentry.  In  many  cases  it  was  only  the  ballistic  coefficient  which  must  be  varied.  This 
sensitivity  of  reentry  time  to  the  ballistic  coefficient  is  given  by  the  following 
approximate  formula  [Ref.  47:p.  39] 
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6*         =  -***  (176) 

reentry  i     i 

« 

where 

tj        =  time  interval  of  motion  prediction  to  the  moment  of  reentry 

k        =  ballistic  coefficient  =  CDA/2m 

Given  a  1  %  change  in  ballistic  coefficient,  the  down  range  travel  is  altered 
by  180  km.  A  2-3%  change  in  ballistic  coefficient  (as  seen  prior  to  the  attempted  control 
maneuver)  alters  the  predicted  impact  point  by  +  500  km.  Figure  48  shows  the 
dependence  of  time,  latitude  and  longitude  on  the  ballistic  coefficient  [Ref.  47:p.  40]. 

The  approximate  estimate  of  the  probable  debris  fallout  range,  considering 
the  errors  in  determining  the  reentry  point,  with  the  vehicle  destruction  occurring  at  75 
km  altitude,  is  equal  to  ±  1000  km  relative  to  the  calculated  impact  point.  This 
corresponds  to  the  sub-satellite  track  crossing  Chili  and  Argentina  from  the  south-west 
to  north-east.  [Ref.  47:p.  40] 
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[Ref.  47] 
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C.      FACTORS  INFLUENCING  REENTRY/IMPACT  DISPERSION 

The  reentry  vehicle's  deceleration  may  be  expressed  in  terms  of  "gravitational 
acceleration  forces"  (g's)  as  follows  [Ref.  96:p.  1] 


Dynamic  pressure        2 ,17-v 

Ballistic  coefficient        W_  . 


Assuming  that  density  and  wind  velocity  are  the  two  atmospheric  characteristics 
which  most  strongly  influence  the  motion  of  a  reentry  vehicle  passing  through  it, 
reference  [96]  states  the  following: 

1 .  Reentry  vehicle  is  non-ablative  and  non-maneuvering  (W  and  A  are  constant) 

2.  Drag  coefficient  is  dependent  upon  altitude  (Mach  dependence) 

3.  Drag  coefficient  is  shape  dependent 

The  displacement  in  impact  point  can  be  considered  to  be  the  sum  of  miss 
contributions  from  the  various  atmospheric  layers,  given  by  the  product  of  influence 
coefficients  times  changes  in  density  or  wind  velocity  for  the  respective  altitude  layers. 
[Ref.  96:p.  2] 

Reference  [97]  studies  the  effects  of  a  non-spherically  symmetric  atmosphere  on 
the  reentry /impact  location  of  a  decaying  satellite  orbit.  This  investigation  deals  with  two 
specific  cases,  a  90  degree  and  30  degree  inclined,  low  Earth  orbit  and  shows  the  effect 
of  the  atmosphere's  diurnal  bulge  on  the  impact  location.  It  is  observed  that  a  decaying 
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object  is  "more  likely"  to  reenter  and  impact  the  Earth  at  certain  latitudes  than  at  others. 
[Ref.  97:p.  iii] 

The  principal  effects  contributing  to  this  impact  distribution,  in  addition  to  the 
diurnal  bulge,  are  the  geometric  oblateness  of  the  atmosphere  and  the  gravitational 
oblateness  of  the  Earth.  Also  identified  as  important  parameters,  are  orbital  inclination 
and  ballistic  coefficient. 

The  results  of  this  investigation  are  as  follows:  [Ref.  97:p.  43] 

1.  Polar  orbit  (90°  inclination)  reentry  is  most  likely  in  the  equatorial  latitudes. 

2.  Latitudes  of  maximum  or  minimum  likelihood  for  impact  are  not  significantly 
affected  by  changes  in  ballistic  coefficient. 

3.  Large,  "balloon  type"  (spherical)  satellites  react  to  the  diurnal  bulge  by  causing 
impact  on  the  opposite  side  of  the  Earth. 

4.  Smaller,  heavier  satellites  show  little  influence  of  the  diurnal  bulge,  but  they 
exhibit  a  greater  variation  in  impact  probability  with  latitude. 

5.  The  impact  distribution  for  a  given  spherical  satellite  "flattens"  as  the  inclination 
is  moved  from  90°  to  30°.  This  is  attributed  to  the  fact  that  the  extremes  in 
Earth  radius  differ  by  only  5  km  under  the  30  °  inclined  orbit  as  opposed  to  21 
km  under  the  polar  orbit. 

Reference  [98],  a  survey  of  satellite  lifetime  and  orbit  decay  prediction,  conducted 

in  1980,  states: 

...gravity  perturbations  cannot,  by  themselves,  lead  to  orbital  decay  since  they  are 
conservative.  But,  they  can  induce  oscillations  in  the  orbit... drag  is  proportional 
to  atmospheric  density... gravity  perturbations  coupled  with  drag  are  more 
significant  than  when  only  drag  is  modeled.  [Ref.  98:p.  3] 
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Also  noted  is  the  observation  that  non-symmetric  mass  distribution  and  shape  may 
amplify  the  torques  created  by:  [Ref.  98:p.  4] 

1.  Gravity  gradients 

2.  Aerodynamic  forces 

3.  Solar  pressure  forces 

An  interesting  note,  also,  is  that  any  mass  loss  increases  the  magnitude  of  the  drag 
deceleration  since  [Ref.  98:p.  5] 

«*»  ~  ~  (178) 

drag       m 

A  final  noteworthy  observation  from  reference  [98]  is  that  any  attempt  to  model 
satellite  attitude,  shape,  lift,  mass  distribution,  ablation  and  breakup  would  largely  exceed 
the  level  of  sophistication  currently  available  (1980)  for  atmospheric  density  models, 
therefore,  these  non-atmospheric  factors  are  neglected  and  treated  as  higher-order  effects. 
[Ref.  98:p.  5] 

Reference  [99]  was  initiated  after  the  reentry  of  part  of  Sputnik  IV  over  Wisconsin 

in  September  1962.  [Ref.  99:p.  2]    This  study  of  factors  influencing  the  prediction  of 

orbit  decay  and  impact  states: 

...final  reentry  occurs  shortly  after  the  point  of  minimum  altitude  (not  coincident 
with  perigee,  because  of  the  Earth's  oblateness)  [Ref.  99:p.  29] 

It  is  possible,  knowing  only  the  initial  inclination  angle  and  reentry  vehicle  orbital 

path,  to  determine  the  points  at  which  the  minimum  altitude  points  will  occur. 
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Therefore,  knowing  where  these  points  occur,  it  may  be  possible  to  predict  reentry  within 
several  tenths  of  an  orbit  revolution,  as  opposed  to  ambiguities  of  one  complete 
revolution  as  noted  in  reference  [95]  (LIFETIME  v3.0),  since  reentry  tends  to  follow  the 
point  of  minimum  altitude.  [Ref.  99:p.  33] 

Without  such  precise  knowledge  of  the  body's  ballistic  coefficient,  the 
characteristics  of  the  atmosphere,  the  orbital  elements  and  their  variation,  only 
"probability  estimates"  can  be  made  of  impact  latitude,  based  simply  on  the  relative  time 
spent  by  the  satellite  in  various  latitude  bands.  This  conclusion  yields  the  following 
relationships:  [Ref.  99 :p.  33] 


1.  If  the  initial  inclination  equals  50°  -  three  times  more  likely  to  impact  between 
40-50°  band  than  0-10°  latitude  band 

2.  If  the  initial  inclination  equals  49°  -  two  times  more  likely  to  impact  the 
continental  U.S.  than  a  65°  inclined  orbit 


A  final  comment  about  factors  influencing  reentry  and  impact  prediction  from  this 

investigation  is  that: 

In  spite  of  the  present  reentry  rate  of  about  one  object  a  day  (most  of  them  burning 
up,  and  2  out  of  3  of  Russian  origin),  the  threat  imposed  on  the  Earth  population 
from  direct  hits  by  debris  from  uncontrolled  reentries  of  space  objects  can 
generally  be  regarded  as  minor  when  compared  with  other  risks  of  civilization 
(traffic  accidents,  etc.).  An  accumulation  of  worst  case  assumptions  could  lead  to 
one  casualty  every  five  years  for  a  densely  populated  area  of  the  size  of  the  Federal 
Republic  of  Germany  in  case  of  no  prior  warning.  In  case  an  early  warning  could 
be  issued,  such  casualties  could  most  likely  be  avoided  completely.  [Ref.  100:p. 
49] 

Figure  49  shows  the  magnitude  of  the  low  Earth  orbit  problem  purely  as  a  function 

of  the  number  of  objects  in  orbit  [Ref.  100:p.  51].   Figure  50  shows  the  limiting  factor 
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Figure  49:  Near-Earth  View  Of  All  Cataloged  Space  Objects  (1987) 
[Ref.  100] 
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Figure  50:  U.S.  Radar  And  Electro-Optical  Assets 
[Ref.  100] 
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of  geographic  location  of  observation  sites  from  which  data  may  be  obtained  for 
reentry /impact  predictions  [Ref.  100:pp.  5-6]. 

The  final  factor  in  reentry/impact  prediction  to  be  discussed  here  is  the  rotation  of 
the  orbital  plane  during  the  final  phase  of  reentry,  as  presented  in  reference  [101].  This 
study  states  that  planar  motion  is  representative  of,  or  an  adequate  model  for,  reentry 
deceleration  and  heating  rate  studies.  It  is,  however,  inadequate  for  the  study  of  ballistic 
reentry  below  altitudes  of  60  km  and  ballistic  coefficients  in  the  ranges  0.001  m2/kg  < 
B  <  0. 1  m2/kg.  For  these  values  of  B,  the  orbital  plane  begins  to  rotate  at  altitudes  of 
30  km  and  60  km  respectively.  [Ref.  101  :p.  1215]  This  investigation  deals  exclusively 
with  the  final  descent  from  120  km,  and  defines  this  as  the  "reentry  phase"  [Ref.  101  :p. 
1216]. 

As  previously  mentioned,  maximum  heating  rates  and  deceleration  occur  at  altitudes 
high  enough,  such  that,  rotation  of  the  orbital  plane  is  negligible  in  the  calculation  of 
these  parameters.  However,  in  the  calculation  of  the  total  range  traveled  to  Earth 
impact,  these  lower  altitude  effects  and  orbital  plane  rotations  become  significant.  [Ref. 
101:p.  1217] 

Several  initial  assumptions  are  made:  [Ref.  101  :p.  1217] 

1.  0.001  m2/kg  <  CDA/m  <  0.1  m2/kg 

2.  0  <  ^o  ^  45°  ,  where  ^0  is  the  initial  reentry  trajectory  angle 

3.  v0  <  15  km/sec 
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As  the  vehicle  descends,  the  ambient  air  becomes  denser  and  the  rotating 
atmosphere  begins  to  force  the  vehicle's  motion  towards  a  direction  parallel  to  that  of  the 
atmosphere's  velocity.  This  effectively  rotates  the  vehicle's  orbital  plane  towards  the 
direction  of  the  atmospheric  velocity  vector.  Figure  51  shows  this  rotational  effect  on 
the  orbital  plane  [Ref.  101  :p.  1219]. 


Figure  51:  Orbital  Plane  Rotation  Due  To  A  Rotating  Atmosphere 
[Ref.  101] 
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V.   STOCHASTIC  AND  STATISTICAL  PREDICTION  METHODS 

This  chapter  will  describe  three  investigations  that  used  stochastic  and  statistical 
techniques  and  methods  to  predict  the  time  and  impact  location  of  decay-induced 
reentry.  Two  investigations  have  direct  applicability  to  the  current  reentry  model  used 
by  the  USSPACECOM  and  the  third  investigation  is  a  Monte  Carlo  simulation  analysis 
that  was  used  to  predict  the  footprint  dispersion  area  of  Skylab. 

A.      EXTENSIONS  OF  THE  PHYSICAL  MODELING  REENTRY  THEORY 

In  Chapters  II  and  IV,  the  current  reentry  theory  used  by  the  USSPACECOM  to 
predict  impact  location  and  time  was  presented.  The  theory's  accuracy  is  limited  by  the 
ability  to  observe  the  reentry  process  and  the  inherent  deficiencies  in  the  model  used  to 
represent  the  physical  reentry  process.  This  section  will  examine  a  Doctoral  dissertation 
and  a  Master's  thesis  from  the  Air  Force  Institute  of  Technology  (AFIT)  which  are 
proposed  extensions  or  modifications  of  the  current  model.  [Refs.  48,102]  These 
investigations  utilize  stochastic  and  statistical  methods  and  techniques  to  improve  the 
predicted  impact  location  and  time. 

1.      Estimation  of  Reentry  Trajectories 

Reference  [102]  developed  a  linearized  differential  corrector  technique  as  an 
extension  of  the  existing  orbital  estimation  technique  into  the  reentry  region  to  estimate 
the  reentry  trajectory.      Reentry  observations  from  a  space  based  sensor  capable  of 
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providing  infrared  angular  trajectory  observations  at  fixed,  discrete  time  intervals  and 

large  uncertainties  in  true  reentry  dynamics,  were  the  primary  engineering  considerations 

of  this  problem  [Ref.  102:p.  1].    The  author  states  the  following  about  the  current 

deterministic  model: 

...the  basic  differential  corrector  with  deterministic  dynamics  is  inadequate  for 
processing  true  reentry  data.  The  deterministic  dynamics  and  infinite  memory  of 
the  basic  formulation,  cause  the  estimator  to  yield  significantly  biased  solutions 
relative  to  the  standard  deviations  of  the  estimator-computed  covariance  matrix. 
These  occur  when  processing  data  reflecting  the  dynamic  variations  anticipated  in 
true  reentry  trajectories.  The  use  of  the  estimator  in  this  form  would  not  provide 
an  accurate  estimate  of  Earth  impacts  for  satellite  debris.  [Ref.  102:p.  19] 

Since  the  uncertain  dynamics  of  the  reentry  process  pose  many  difficulties  for  the 

existing  deterministic  model,  a  significant  portion  of  the  research  in  the  dissertation  was 

devoted  to  identifying  limits  of  various  estimation  techniques  that  were  considered  for 

this  problem.  [Ref.  102:p.  1] 

An  adaptively  determined,  ad  hoc,  scalar  finite  memory  or  fading  memory 
parameter  approach  was  selected.  The  motivation  of  this  approach  is  based  on  the  fact 
that  the  existing  reentry  theory  attempts  to  group  all  unmodeled  parameters  as  a  constant 
with  the  ballistic  coefficient  (CDA/m)  over  some  short  trajectory  span.  [Ref.  102:p.  11] 

In  order  to  extend  the  existing  orbit  determination  technique  into  the  satellite 
reentry  problem,  reference  [102]  defined  a  reentry  dynamics  model  and  developed  the 
weighted  least-squares  differential  corrector  structure  currently  used  in  the  deterministic 
model.  The  estimator  update  and  propagation  equations  of  the  differential  corrector 
structure  for  an  infinite  memory  formulation  are  summarized  as  follows: 
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Propagations  between  epochs: 

State:   integrate  with  initial  conditions  x^O,  from  the  epoch  at  t^ 

m  =  mm  (179) 

to  obtain  xm(tm  +  ,)  at  epoch  tm  +  l 
Covariance: 

sm,ljn  -  *(»„.„«»)  s„,„  *(»„.„',/  (180) 

Update  at  next  epoch: 
State: 

/ 

*m*l(fm+d   =  Xm<fm*l)   +  E  8j#-l)|  (181) 

i-1 

Covariance: 

$..,,.,♦,  =  (i..,,.-'  ♦  V  VV  (182) 

Recall  that  1  is  the  number  of  iterations  required  to  satisfy  convergence  and  T(n)  T 
R(n)"1  T(n)  is  evaluated  from  the  reference  trajectory  on  the  final  iteration,  1.  [Ref. 
102:pp.  33-34] 

Equation  (179)  is  a  nonlinear  dynamics  model,  where  f(x(t),  t)  is  a  deterministic  function 

of  the  state  variables  and  is  a  continuous  function  of  time.    The  variable,  5x(t),  from 

equation  (181)  is  a  "seeked"  correction  that  will  minimize  a  weighted  quadratic  cost 

function  of  the  observation  residuals.    Iterations  of  the  differential  corrector  continues 

until  the  observation  residuals  convergence  criteria  is  satisfied.    Equation  (180)  is  the 

initial  state  covariance  matrix  that  remains  constant  until  the  iterative  process  has 

converged.    The  quantity  <£(tm+i>tm)  from  equation  (180),  is  the  state  transition  matrix. 

The  final  iteration  in  the  process  yields  equations  (181)  and  (182).    R(n)and  T(n)  from 

equation  (182)  are  respectively  the  matrices  of  the  observation  noise  covariances  and  the 

observation  index  being  processed.  [Ref.  102:pp.  20-29]  As  previously  mentioned: 
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Application  of  this  batch  processing  algorithm  to  reentry  estimation  has  often 
resulted  in  poor  estimator  performance.  This  is  largely  due  to  the  more  significant 
non-linear  dynamics  of  reentry  and  the  use  of  a  deterministic  and  simplistic 
dynamics  model  ...  in  an  uncertain  dynamics  region.  [Ref.  102:p.  31] 

In  the  development  of  the  reentry  dynamics  model,  the  author  chose  an  eight 

dimensional  state  vector  defined  as   [Ref.  102 :pp.  41-42] 


xl  =  x 


x2  =  x 


x(t)  = 


x3  =  y 

X4  =  > 


x5  =  z 


x6  =  z 

*7    =  BP0 
*i   =   <? 


(183) 


where 


x*y>z  =  velocity  components 

x,y,z,  =  position  components  (cartesian) 

B  =  ballistic  coefficient  (CDA/m) 

p0  =  sea  level  air  density  from  equation  (10) 

Q  =  density  scale  height   (RT0  /  Mog0*)  from  equation  (10) 

Mo  =  molecular  weight  of  air 

g0*  =  acceleration  of  gravity  at  Earth  surface 

T0  =  atmospheric  temperature  at  sea  level 

R  =  gas  constant 


(184) 
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H*      =  altitude  (g/goHo)  from  equation  (10) 

g        =  local  gravitational  acceleration  from  the  geopotential  model  x,y,z 
coordinates 

go  =  reference  geoid  level  gravitational  acceleration 

H0  =  geocentric  altitude  (R-RJ 

R  =  local  radius  position  relative  to  the  Earth  center  (x2  +  y2  +  z2)1/2 

R,  =  Rq(1  -f)[l  -(2f-f5)cos25]1/2 

f        =  flattening  factor  of  the  reference  geoid  whose  elliptical  shape  is  consistent 
with  J2 

5        =  local  latitude  (  cos^  (x2  +  y2)m  I  Ro]  ) 

Ro      =  radius  of  the  reference  geoid  at  the  equator 

The  exponential  atmospheric  density  model,  equation  (10),  was  used  in  the 

reentry  estimator  because  of:  [Ref.  102:p.  37] 

1.  The  reduced  mathematical  complexity. 

2.  The  availability  of  continuous  valued  density  and  partial  derivatives  of  the 
density  along  the  reentry  trajectory. 

The  eighth  state  variable,  Q,  was  chosen  since  it  is  slowly  changing  at  the 
reentry  altitudes  (  <  100  km).  For  this  application,  the  initial  value,  Q  =  7.0031  km 
and  its  covariance  were  derived  using  a  least  squares  fit  to  the  base  density  values  of  the 
altitude  layers  from  the  U.S.  Standard  Atmosphere  1962.  Also,  an  Earth  Centered 
Inertial  (ECI)  coordinate  frame  was  chosen  to  minimize  the  complexities  of  the 
observation  noise  covariance  matrix,  R(n).  [Ref.  102:p.40] 
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Both  aerodynamic  and  geopotential  terms  within  the  deterministic  function  of 
state  variables,  f(x(t) ,  t),  are  included  in  the  estimator  dynamics  model.  The  dynamics 
equation  is  given  in  the  form  [Ref.  102:pp.  42-44] 


x  = 


*1 

=  *2 

*2 

=4 

+  4 

*3 

=  *4 

*4 

=4 

+4 

*3 

=  *6 

*6 

■A 

+4 

*1 

=  0 

*8 

=  0 

(185) 


The    aerodynamic  acceleration  x,y,z  components  from  equation  (185)  are 
derived  from  equation  (6) 


4  =  -\BpVA(x  +  ay) 


(186) 


4  =  -^BpK^^-wx) 


(187) 


4  *  4Bp  ^  * 


(188) 


The  geopotential  acceleration  terms  from  equation  (185)  were  derived  from 
the  Smithsonian  Astrophysical  Observatory  SAO-III  Earth  Model  using  the  model 
parameters  and  zonal  and  tessera!  coefficients 
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/_    =  coscc.KK,  -  sina  W„ 
f    =  sina  Wx  +  cosaaWv 

By  m         *  s  j 

4 '  ^ 

where 

a,       =  Right  Ascension  (RA)  of  the  Greenwich  Meridian 
The  gradient  terms  from  equations  (189)  through  (191)  are  defined  as 


(189) 
(190) 
(191) 


^  =  EE 


(        BU?  dVn 

c„  — -  +  s. 


n=Q    m=0 


m\ 


dc 


&  J 


(192) 


m\ 


A   A   I         31/"  dV. 

n-o  m=o  V       3y  3y  / 


(193) 


^  =  E  E 

n=0     m=0 


C._  — -  +  5. 


dz 


dz) 
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where 


Cnm     =  zonal  harmonic  coefficient 


Sno,     =  tesseral  harmonic  coefficient 


_       GM  Rq 
UT  =  —   -5-  PBw(6/)  cos(wX) 


tfB+1 


(195) 
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,  * 


V;  =  ^A  P;(l)  sin(mA)  (196) 

■  nn+l  ■ 


G       =  universal  gravitational  constant 

M       =  mass  of  Earth 

Ro      =  mean  equatorial  Earth  radius 

R       =  distance  from  Earth  center 

Pnm     =  associated  Legendre  polynomial  of  degree  n  and  order  m 

5'       =  z/R   (z  =  coordinate) 

X        =  longitude 

As  previously  mentioned  in  Chapter  II,  in  the  region  below  120  km,  the  first-order  zonal 

harmonic,  J2,  may  attain  a  magnitude  approaching  that  of  the  atmospheric  drag.   Based 

on  this  fact,  in  the  development  of  the  actual  estimator  and  simulation,  only  the  central 

gravity,  C^,  and  Earth  oblateness  (J2),  C20 ,  terms  were  used  in  the  reentry  altitude 

regions. 

Numerical  simulations  using  simulated  reentry  data  derived  from  a  realistic 

"truth  model"  were  conducted  on  the  basic  estimator  structure  and  dynamics  model  to 

quantify  its  performance. 

These  analyses  examined  the  effects  of  mismatch  between  the  truth  model  and  the 
estimator  model  dynamics,  accuracy  variations  on  the  angular  observations, 
multiple  orbital  observation  locations  and  variations  of  the  geometric  relationships 
between  the  observing  satellites  and  the  reentry  trajectory.  [Ref.  102:p.  14] 
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A  series  of  Monte  Carlo  analyses  were  conducted  which  showed  the  model  dynamics  as 
the  most  important  factor  impacting  the  estimator  performance.  This  was  accomplished 
by: 


1.  Considering  the  velocity  estimates  and  mean  bias  of  the  trajectory  position  vs 
the  magnitude  of  the  standard  deviations  from  the  estimator  computed  state 
covariance  matrix. 

2  Comparing  the  magnitudes  of  the  standard  deviations  from  the  estimator 
computed  state  covariance  matrix  to  those  derived  in  the  Monte  Carlo  samples. 


In  the  search  to  find  possible  solutions  to  the  problems  associated  with  the 
estimator  performance,  reference  [102]  reviewed  and  discussed  the  limitations  of  several 
model  compensation  methods  that  are  relative  to  the  reentry  application,  including:  [Ref. 
102:p.  15] 

1.  Adding  a  pseudo-noise  compensation  to  the  model  dynamics 

2.  Adaptive  estimation  methods 

3.  State  covariance  deweighting  techniques 

The  author  developed  a  fading  memory  differential  model  compensation  method  based 

on  the  fact: 

. .  .the  previous  numerical  results  indicated,  the  fundamental  limitation  of  the  infinite 
memory  estimator  formulation  with  deterministic  dynamics  is  the  biased  estimator 
solutions  which  occur  when  processing  true  reentry  dynamics.  With  i)  exact 
dynamics  ii)  an  upper  limit  on  the  time  span  of  valid  linearization,  and  iii)  a 
lower  limit  on  observation  accuracy  (greater  than  or  equal  to  10"5  radians), 
acceptable  estimator  performance  is  available  in  terms  of  bias  and  RSS/(ONE 
SIGMA)  Ratio.  [Ref.  102:p.  112] 
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Since  the  variations  in  global  atmospheric  density  and  the  dynamic  changes 
of  the  reentry  body  pose  such  a  virtually  intractable  problem  within  a  deterministic  set, 
improving  the  linearized  estimator  performance  is  based  on  the  impact  of  uncertainties 
within  the  ballistic  coefficient  and  the  atmospheric  density  [Ref.  102:pp.  112-113]. 
Improvement  of  the  estimator  was  achieved  by  an  adaptive  determination  of  an  ad  hoc 
scalar  multiplier,  7.  A  finite  memory  on  the  processing  of  earlier  observations  is 
implemented  by  multiplying  the  parameter,  7,  to  the  terms  of  the  state  covariance  matrix 
prior  to  an  observation  update.  At  each  update  along  the  trajectory,  the  size  of  the 
change  in  state  variable,  5x;  is  examined.  Each  5x{  magnitude  is  then  compared  to  the 
magnitude  of  the  standard  deviation  of  their  respective  terms  within  a  "deweighted" 
covariance  matrix  which  is  computed  by  the  estimation  algorithm.  As  with  the  infinite 
memory  estimator  model,  a  series  of  Monte  Carlo  analyses  were  conducted  to  assess  the 
model's  ability  to  estimate  an  anticipated  true  reentry  dynamics  trajectory.  [Ref.  102:pp. 
15-16] 

Finally,  reference  [102]  demonstrated  the  ability  of  the  fading  memory 

method  to  provide  a  tangent  plane  projection  of  Earth  impact  locations  as  shown  in 

Figure  52  [Ref.  102:p.  192],  by  using  bias  magnitudes  within  the  standard  deviations  of 

the  deweighted  state  covariance  matrix. 

The  standard  deviation  of  the  position  error  from  the  deweighted  state  covariance 
matrix  provides  a  good  definition  of  the  uncertainties  in  the  estimated  Earth  impact 
locations,  thus  it  can  be  used  to  define  a  search  area  for  recovery  of  satellite 
debris.  These  results  illustrate  the  viability  of  the  method  to  estimate  decayed 
satellite  impact  locations  and  uncertainties  significantly  improved  over  existing 
astrodynamic  applications.  [Ref.  102:p.  16] 
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Figure  52:  Tangent  Plane  Projection  Of  Earth  Impact  Location 
[Ref.  102] 
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The  primary  conclusions  and  recommendations  from  the  dissertation  are: 
[Ref.  102:pp.  198-202] 


1.  Dynamics  uncertainties  of  the  general  satellite  decay  trajectories  significantly 
affect  the  estimator  performance  as  shown  by  the  Monte  Carlo  simulation 
analyses. 

2.  Angular  observation  accuracies  with  standard  deviations  less  than  10"5  radians 
induce  significant  error  in  the  state  estimate  vs  the  standard  deviations  of  the 
state  covariance  matrix  in  deterministic  dynamics  models. 

3.  A  recursive  formulation  of  the  estimator  is  recommended  that  uses  a  short  time 
span  between  the  epoch  or  trajectory  update  point  and  the  observation(s)  being 
processed  due  to  the  anticipated  dynamics  uncertainties  of  reentry. 

5.  Multiple  observations  from  more  than  one  the  orbital  source  are  required  to 
improve  the  observability  of  the  reentry  and  to  provide  higher  data  content  over 
similar  time  spans  in-order  to  achieve  acceptable  estimator  performance  in  terms 
of  bias  and  RSS/(ONE  SIGMA)  ratio. 

6.  An  eight  dimensional  state  vector  provided  superior  estimation  performance  as 
compared  to  the  seven  dimensional  vector  used  in  the  current  deterministic 
models.  Performance  improvement  is  achieved  through  simpler  mathematics  in 
the  dynamics  model  and  continuous  valued  partial  derivatives  of  the  dynamics 
over  the  trajectory  space  for  a  Taylor's  series  linearization. 

7.  The  adaptive, ad  hoc  scalar  fading  memory  parameter  is  easily  incorporated  into 
the  basic  estimator  structure. 

8.  A  Monte  Carlo  derived  impact  covariance  for  the  final  propagation  phase  to 
impact  preserves  the  integrity  solution  statistics  over  the  non-observable  final 
portion  of  the  trajectory.  Impact  location  uncertainties  are  on  the  order  of  one 
to  two  magnitudes  smaller  than  those  available  from  current  operational 
techniques. 

9.  Further  investigations  which  vary  the  observational  data  rate  and  the  time 
varying  character  of  true  reentry  dynamics  are  needed  to  examine  the  estimator 
performance  extensions.  Specifically,  more  accurate  observations  (much  less 
than  10'5  radians)  with  higher  data  rates  and  frequency  variations  as  well  as 
alternative  measurements,  such  as  range  or  range  rate. 
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10.  Further  analysis  of  applying  the  fading  memory  method  to  very  high  altitudes 
and  shallow  reentry  angles  is  recommended  due  to  the  fact  that  violent  dynamic 
changes  under  these  conditions  may  result  in  divergent  estimator  performance. 

11.  Application  of  this  estimation  technique  to  a  wide  class  of  reentry  trajectories 
could  provide  a  large  empirical  data  base  to  improve  the  estimator  dynamic 
model.  Dynamic  model  pseudo-noise  compensation,  statistical  linearization,  or 
higher  order  filters  investigations  should  be  pursed. 

2.       Analysis  of  Tracking  and  Impact  Prediction  (TIP) 

Reference  [48]  analyzed  the  accuracy  of  early  TIP  processing  (Chapter  II, 
p.58)  conducted  by  the  USSPACECOM  on  180  objects  that  decayed  during  the  years 
1987  through  1990.  As  part  of  the  analysis,  early  TIPs  (7  day  to  3  hour  predictions) 
were  compared  to  the  final  TIP.  The  time  error  for  each  TIP  run  was  calculated  and 
compared  to  the  ±20%  accuracy  level  claimed  by  the  SSC  at  Cheyenne  Mountain  AFB. 
Results  from  the  comparison  study  indicate: 


1.  The  decay  prediction  accuracy  is  usually,  but  not  always  within  the  claimed 
accuracy  level  as  shown  in  Figure  53  [Ref.  48  :p.  37]. 

2.  The  existence  of  a  positive  bias  which  indicates  early  TIPs  are  routinely  late 
relative  to  the  final  TIP. 


The  author  developed  six  multiple  linear  regression  models  that  could  be  incorporated 
into  the  TIP  decay  procedures  in  an  attempt  to  model  out  some  of  the  positive  bias  found 
in  TIP  decay  prediction  data.  [Ref.  48:p.  viii] 
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Figure  53:  Decay  Predicted  Accuracy  (By  Year) 
[Ref.  48] 
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The  data  used  in  the  investigation  was  collected  from:  [Ref.  48:pp.  15-17] 


1 .  TIP  Required  Item  Checklist~a  manually  kept  chronological  account  of  each 
TIP  processing. 

2.  Decay  History ~a  computer  generated  log  of  each  SP  update  that  includes  the 
run  time,  time  of  last  observation,  epoch  time,  epoch  revolution  number,  B-term 
(ballistic  coefficient),  period  and  decay  prediction. 

3.  Final  Element  Set~a  listing  of  the  final  orbit  parameters.  The  eccentricity  and 
mean  motion  data  from  the  listing  were  used  in  the  investigation. 


Only  those  objects  that  received  the  entire  TIP  update  cycle  (7  day  through 

final  run)  were  used  to  accurately  analyze  the  effects  of  each  successive  update  and 

prediction  [Ref.  48:p.  17].    In  order  to  assess  the  accuracy  level  claimed  by  the  SSC, 

reference  [48]  used  the  data  from  the  final  decay  prediction  as  the  control  by  which  to 

compare  the  earlier  predictions. 

The  final  prediction  was  chosen  as  the  control  because  it  uses  observations  which 
are  closest  to  the  actual  impact  point  and  is  considered  to  be  the  most  accurate 
prediction  available.  In  order  to  further  justify  the  use  of  the  final  prediction  as 
the  control,  a  statistical  analysis  was  also  performed  to  directly  compare  the  few 
sighted  reentry  points,  called  Vis  Obs,  with  the  final  prediction  made  by  the  Space 
Surveillance  Center.  [Ref.  48:p.  18] 

Of  the  180  objects  studied,  93  were  Vis  Obs.  Figure  54  shows  the  mean  time 

error  of  the  final  run  time  vs  the  Vis  Obs  [Ref.  48:p.25].  The  size  of  the  final  time  error 

standard  deviation  decreases  as  shown  in  Figure  55  [Ref.  48:p.  26].  This  decrease  may 

be  correlated  to  the  level  of  solar  activity  during  the  1987-1990  time  period.    Solar 

activity  levels  began  to  dramatically  increase  in  1987  and  continued  to  increase  through 

the  solar  maximum  (March  1990).   The  rate  of  change  increase  of  the  sunspot  activity 
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Figure  54:  Final  Run  vs  Vis  Obs  Mean  Time  Error 
[Ref.  48] 
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Figure  55:  Final  Time  Error  Standard  Deviation 
[Ref.  48] 
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and  solar  flux  during  1987-1988  was  greater  than  in  period  from  1989-1990.  According 
to  the  author,  the  trend  shown  in  Figure  55  may  be  related  to  the  lesser  rate  of  change 
during  1989-1990.  [Ref.  48:pp.  26-28] 

The  time  error  for  each  separate  TIP  run  was  calculated  as  the  difference 
between  the  predicted  decay  time  for  that  run  and  the  final  run  [Ref.  48:p.  18].  Figure 
56  shows  the  graphic  calculation  results  for  the  mean  time  error  and  time  error  standard 
deviation  [Ref.  48:p.  29]. 

In  addition  to  the  time  error  calculations,  the  location  error  was  calculated 
by  taking  the  difference  between  the  predicted  location  point  for  each  run  and  the  final 
predicted  point.  The  method  used  the  mean  motion  (n)  from  the  Final  Element  Set  to 
accurately  determine  the  velocity  for  each  TIP  object.  Since  mean  motion  was 
unavailable  for  the  early  TIP  runs,  an  approximation  was  made  using  the  final  mean 
motion  value  for  each  TIP  run.  This  introduced  an  error  into  the  calculation  of  the 
location  error.  However,  because  the  location  errors  are  very  large  (thousands  of  km), 
the  error  introduced  by  using  the  final  mean  motion  instead  of  the  actual  mean  motion 
for  that  particular  TIP  run  was  considered  insignificant.  The  mean  motion  value  was 
used  to  calculate  the  semi-major  axis  distance  from  the  equation 
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3  (197) 


The  velocity  was  calculated  using  the  semi-major  axis  from  the  equation 


207 


MEAN  TIME  ERROR 
(1987-1990) 


*»/ 


40> 


40U3 


T-OAY       «SAY   '  MAT      1<CAT       »H»        «-WL        y»k 
TIP  RUN 


TIME  ERROR  (1987-1990) 
(Standard  Deviation) 


7-CAY    «AY     WAY     1-OAY     U+iR      (Mffi       Mtt 
TV  RUN 


Figure  56:  Mean  Time  Error  And  Standard  Deviation  (1987-1990) 
[Ref.  48] 
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The  object's  velocity  was  then  multiplied  by  the  previously  calculated  time  error  in  order 
to  determine  the  location  error  (difference  between  the  predicted  location  point  for  each 
run  and  the  final  predicted  point).  [Ref.  48:pp.  19-22]  Figure  57  shows  the  graphic 
calculation  results  for  the  mean  location  error  and  location  error  standard  deviation  [Ref. 
48:pp.  32-33]. 

Multiple  linear  regression  describes  the  relationship  between  several 
independent  variables  and  a  dependent  variable.  The  motivation  for  developing  the 
multiple  linear  regression  model  to  eliminate  the  bias  found  in  the  mean  time  error  was 
based  on  the  second  objective  of  the  investigation.  This  objective  was  to  determine  if 
it  was  advantageous  to  initiate  an  OPREP-3  report  (used  to  notify  higher  authority  of  a 
potential  reentry  within  100  nm  of  the  former  Soviet  Union  border)  earlier  than  the  6 
hour  mark.   The  first  step  in  the  process  was  to  determine  if  the  model  in  the  form 

£fc)    =    Po   +    Pl'l    +    Plh    +    P3f3    +    PA    +    P5'5    +    P6'6   +    P7'7  (199) 

where 

E(tf)  =  expected  value  of  the  final  decay  prediction  time 

tj  =  early  TIP  decay  predictions 

tf  =  final  decay  prediction  time 

ft  =  y-intercept  and  coefficients  to  be  determined 
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Figure  57:  Mean  Location  Error  And  Standard  Deviation  (1987-1990) 
[Ref.  48] 
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could  be  found  where  the  early  TIP  data  could  be  used  to  approximate  the  final  decay 

prediction  time.  [Ref.  48:pp.  22-23]    The  author  states  the  following  about  equation 

(199): 

...multiple  linear  regression  was  then  used  to  first  determine  if  a  model  in  the  form 
shown... could  be  found  to  predict  the  final  decay  time.  The  results  were  an  R- 
squared  value  of  1.0000  and  a  p-value  of  .0001.  This  means  that  at  a  significance 
level  of  .05  there  exists  a  perfect  linear  relation  between  some  of  the  independent 
variable  and  the  dependent  variable  where  at  least  two  of  the  /S  terms  are  not  zero. 
The  variance  inflation  values  were  all  extremely  large,  indicating  the  independent 
variables  were  all  highly  correlated  and  that  a  great  deal  of  redundancy  exists  in 
the  data.  [Ref.  48:pp.  37-38] 

Based  on  the  above  results,  six  separate  linear  models  were  developed  in  an  effort  to 

approximate  the  final  TIP  with  greater  accuracy  than  the  current  process.     Specific 

characteristics  used  in  the  modeling  include  the  following  :  [Ref.48:p.  38] 

1.  The  first  model  uses  only  7-day  TIP  data  to  calculate  the  expected  value  of  the 
final  decay  prediction  time,  E(tj). 

2.  One  additional  decay  prediction  data  point  is  incorporated  in  each  subsequent 
model. 

3.  All  180  TIP  objects  were  used  in  the  six  models  to  calculate  E(V). 

The  six  models  are  defined  as 
EJtjj  =  -0.116442  +  0.999064ft)  (200) 

EJtj)  =  -.0155478  +  0.49007ft)  +  0.951 197ft.)  (201) 

EjtJj  =  -.0082444  +  0.41692ft)  +  0.050001ft)  +  0.0908343ft)  (202) 
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EJtJj  =  -0.053030  +  0.007471(f1)  -  0.014395(f2)  +  0.196518(f3)  (2Q3) 

+  0.810492(f4) 

EJtjj  =  -0.022322  -  0.004222^)  -  0.005502(f2)  +  0.046458(f3)  (2Q4) 

-0.049000(f4)  +  1.0123  ll(r5) 

EJtjj  =  -0.008370  +  0.000553^)  -  0.002944(f2)  +  0.014541(f3)  (2Q 

-0.009096(r4)  -  0.183413(r5)  +  1.180378(f6) 

where 

tj        =  7-day  prediction  time 

tj        =  4-day  prediction  time 

t3        =  2-day  prediction  time 

t4        =  1-day  prediction  time 

t5        =  12-hour  prediction  time 

t$        =  6-hour  prediction  time 

Figure  58  shows  the  mean  approximate  error  for  the  six  regression  models 
[Ref.  48:p.  40].  By  comparing  Figure  58  with  the  mean  time  error  in  Figure  56,  the 
regression  models  show  a  better  mean  approximation  error  than  the  TIP  runs.  The 
conclusions  and  recommendations  from  the  thesis  are:  [Ref.  48:p.42] 

1 .  The  decay  predictions  were  much  better  in  general  than  the  reported  ±  20  % . 

2.  The  use  of  linear  models  in  conjunction  with  the  data  generated  by  the  TIP 
processing  would  allow  the  SSC  to  better  predict  final  decay  time  by  elimination 
of  positive  bias  in  the  data. 
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Figure  58:  Regression  Model  Mean  Approximate  Error  (1987-1990) 
[Ref.  48] 
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3.  The  SSC  should  conduct  a  study  of  the  current  Special  Perturbations  model  to 
attempt  to  better  account  for  the  level  of  solar  activity  and  its  effect  on  the 
atmosphere. 


B.   MONTE  CARLO  ANALYSIS  OF  SKYLAB'S  IMPACT  AREA 

Reference  [45]  investigated  the  debris  impact  point  dispersion  area  of  Skylab.  A 
combination  of  Monte  Carlo  statistical  analysis  and  parametric  methods  were  used  to 
determine  the  three-sigma  limits  of  the  debris  footprint.  The  investigation  was  conducted 
by  Martin  Marietta  Aerospace  during  the  design  and  development  of  the  Teleoperator 
Retrieval  System  (TRS)  for  the  Skylab  reboost/deboost  mission.  The  dispersion  analysis 
was  conducted  to  support  the  deboosting  of  Skylab  to  a  safe  oceanic  impact  area  clear 
of  islands  and  routine  shipping  lanes.  [Ref.  45:pp.  1-2] 

The  Monte  Carlo  statistical  analysis  is  an  efficient  and  realistic  approach  that  can 
be  used  to  calculate  the  debris  impact  area  because  of  the  large  number  of  input  variables 
and  nonlinearities  of  the  problem.  Impact  dispersion  area  boundaries  depend  on  the 
following  factors:  [Ref.  45: p.  2] 

1.  Entry  dispersions 

2.  Relative  flight  path  angle 

3.  Debris  ballistic  coefficient 

4.  Breakup  altitude 

5.  Environmental  conditions  (wind  direction/magnitude  and  atmospheric  density) 
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The  entry  dispersions  were  represented  by  a  6  x  6  covariance  matrix  of  velocity 
and  position  errors  at  the  entry  altitude.  For  the  deboost  mission,  the  error  sources  were 
the  uncertainties  associated  with:  [Ref.  45 :p.  2] 

1 .  Accelerometers 

2.  Gyros 

3.  Computer  errors 

4.  Initial  alignment  errors 

5.  Vehicle  performance  dispersions 

The  covariance  matrix  of  state  variables  was  determined  by  conducting  a  trajectory  error 
analysis  using  a  1-sigma  deviation  of  the  errors  and  the  nominal  flight  trajectory.  In  the 
application  of  the  Monte  Carlo  analysis  the  major  considerations  were:  [Ref.  45 :p.  2] 

1.  The  simulation  of  perturbed  trajectories. 

2.  The  modeling  of  the  satellite  breakup. 

3.  The  determination  of  impact  points  for  each  trajectory. 

At  the  entry  altitude,  the  Monte  Carlo  method  required  generation  of  random  state 
vectors  from  the  6  x  6  covariance  matrix.  Because  the  velocity  and  position  errors  are 
correlated,  the  random  error  vector  in  an  uncorrelated  space  was  derived  from  a 
transformation  to  a  principle  axes  system.   This  vector  is  given  by 
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er  -  [nxev  nfy  ....  n6ej 


(206) 


where 

rij       =  a  set  of  random  numbers  drawn  from  a  normal  distribution  of  mean  0  and 
variance  1 

e{       =  the  square  roots  of  eigenvalues  of  the  covariance  matrix 

The  original  coordinate  system  random  state  vectors  were  derived  from  the  equation 


I* 


(207) 


where 


4>        =6x6  eigenvector  matrix  that  transforms  the  principle  axes  system  into 
the  original  coordinate  frame 

Rp      =  perturbed  radius  vector 

Vp      =  perturbed  velocity  vector 

R„      =  nominal  radius  vector 

Vn      =  nominal  velocity  vector 
In  order  to  calculate  the  dispersion  area,  the  Monte  Carlo  analysis  required  a  large 
number  of  simulated  trajectories  and  random  perturbed  state  vectors  at  the  entry  point 
for  each  perturbed  initial  condition.  [Ref.  45:pp.  2-3] 

Additionally,  the  Monte  Carlo  scheme  required  a  fast  and  efficient  computational 
technique  to  calculate  flight  trajectories  from  several  initial  conditions.  Analytical 
solutions,  such  as  Loh's  second-order  theory  (presented  in  Chapter  III)  were  not  suitable 
for  footprint  dispersion  analysis  because  these  solutions  are  valid  for  only  certain  portions 
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of  the  trajectory  and  entry  conditions.    The  trajectory  simulation  used  the  equations  of 
motion  of  a  nonlifting  vehicle  in  the  form 

v  pV(vx  +  <»y)g  (208) 

*      **  2BP 


V   =  o    -   pV(Vy  +  ^g  (209) 

y      *>  2BP 


V    -  a     -    PK(K^  (210) 

2      **  2BP 

where 


V  =  J{VX  +  ay)2  +  (Vy  -  (Six)2  +  (Vf  (relative  velocity)  ^211> 

Sx>gy>gz     =  gravitational  acceleration  components 

Vx,Vy,Vz  =  inertial  velocity  vector  components 

x,y,z        =  radius  vector  components 

o)  =  Earth's  rotation  rate 

BP  =  ballistic  parameter  or  coefficient  (W/CDA) 

Equations  (208)  through  (210)  were  integrated  several  times  from  entry  to  the  impact 
point  for  the  Monte  Carlo  analysis.  The  integration  step  size  used  was  as  large  as 
possible  due  to  a  stability  consideration.  The  system  approached  a  dynamic  equilibrium 
condition  which  resulted  in  a  numerical  instability.  This  was  a  result  of  the  opposing 
gravitational  acceleration  and  aerodynamic  deceleration  during  the  vertical  portion  of  the 
flight.    A  numerical  search  was  conducted  to  determine  the  step  size  as  a  function  of 
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satellite  debris  ballistic  coefficient.  [Ref.  45:pp.  3-4]    Figure  59  shows  the  stable  and 

unstable  regions  used  for  the  numerical  integration  [Ref.  45 :p.  4]. 

Since  the  Monte  Carlo  analysis  required  the  simulation  of  a  large  number  of 

trajectories,  an  accurate  computationally  fast  atmospheric  density  model  was  needed. 

The  authors  chose  the  1962  U.S.  Standard  Atmosphere  that  provided  density  values  from 

0  to  400,000  ft.  [Ref.  45:p.  5] 

Although  both  static  and  dynamic  global  atmospheric  models  accounting  for 
diurnal,  seasonal,  and  latitudinal  variations  are  available,  it  is  found  that  such 
models  are  not  suitable  for  the  Monte  Carlo  dispersion  analysis  of  satellite 
footprints.  A  sensitivity  study  of  density  variation  indicates  only  a  second  order 
effect  on  footprint  dispersion.  Furthermore,  the  results  of  Purcell  and  Barbery 
show  the  downrange  impact  errors  resulting  from  atmospheric  variations  are  slight. 
[Ref.  45:p.  5] 

The  simulation  assumed  one  breakup  altitude  where  all  the  pieces  had  the  same 

initial  position  and  velocity  at  breakup.    By  using  a  bounded  parametric  approach,  the 

smallest  and  largest  debris  ballistic  coefficients  were  used  to  determine  the  largest 

uprange  and  downrange  footprint  dispersions  from  the  nominal  impact  point.    Since  the 

ballistic  coefficient  will  vary  as  it  passes  through  the  atmosphere  due  to  variations  of  CD 

and  Mach  number,  the  program  used  a  generalized  CD  vs  Mach  number  curve  for 

tumbling  pieces  to  update  the  ballistic  coefficient  at  various  altitudes.    This  curve  is 

based  on  the  results  of  a  range  safety  study  of  Titan  launch  vehicle  debris.  [Ref.  45  :p. 

5] 

For  the  footprint  dispersion  study,  it  is  found  that  the  transonic  flow  region,  where 
CD  variation  with  Mach  number  is  significant,  occurs  during  the  vertical  descent 
of  the  debris  with  minimum  effect  on  footprint  dispersion.  However,  the  variation 
of  ballistic  parameters  with  altitude  affects  the  time  of  arrival  of  the  debris  on  the 
ground.  [Ref.  45 :p.  5] 
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Figure  59:  Numerical  Stability  Region 
[Ref.  45] 
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The  latitude  and  longitude  (lat/long)  impact  points  of  the  smallest  and  largest 

Skylab  debris  were  provided  by  the  trajectory  simulation  for  the  various  initial 

conditions  at  the  entry  altitude.    The  Lat/Long  impact  points  were  statistically  treated 

which  determined  a  three-sigma  impact  area  boundary.     The  point  estimate  and 

confidence  level  estimate  of  the  impact  points  for  any  given  confidence  level,  e,  and 

probability,  a,  are  the  fundamental  statistical  quantities  of  interest. 

These  statistical  computations  relate  the  calculated  longitude  and  latitude  bounds 
from  a  finite  sample  to  their  true  values  corresponding  to  an  infinite  sample 
required  for  the  Monte  Carlo  technique.  [Ref.  45 :p.  6] 

By  using  a  normal  distribution  approximation  given  by 

(212) 


N  - 

olK] 

(1 

-  «) 

where  K€ 

is  defined  from  the 

equation 

equation  (212)  provides  the  estimated  number  of  samples  required  to  estimate  a  good 
point.   A  probabilistic  statement  given  by 

Pr(0  *  0B)  =  a  (214) 

Pr®**J  =  a  (215) 

where 

4>       =  impact  point  latitude 
8        =  impact  point  longitude 
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was  obtained  from  the  frequency  histogram  and  cumulative  probability  distribution  of 
Lat/Long  impact  points  generated  from  the  Monte  Carlo  analysis.  For  any  specified 
value  of  a,  the  object  of  the  statistical  analysis  was  the  estimation  of  (/>„  and  6a.  [Ref. 
45  :p.  6] 

1.       Simulation  Results 

During  the  simulation  to  determine  the  debris  footprint  dispersion,  Skylab  was 

deboosted  from  a  170  nm  circular  orbit.    The  nominal  entry  and  corresponding  impact 

point  Lat/Long  were  respectively  (30°S,30°W)  and  (49.86°S,34.92°E).    The  nominal 

entry  point  state  vector  and  the  lower  half  of  a  symmetrical  ECI  frame  covariance  matrix 

are  shown  in  Table  XVI  [Ref.  45  :p.  7].  By  using  a  Gaussian  distribution  for  the  Monte 

Carlo  analysis,  the  initial  state  vectors  were  randomly  selected  from  the  covariance 

matrix.     Specifically,  500  random  entry  states  were  created  by  using  3000  random 

numbers  drawn  from  the  distribution. 

To  ensure  that  these  numbers  truly  represent  a  normal  distribution,  their  mean  and 
variances  were  determined  and  adjusted  to  be  0  and  1  within  stringent  tolerances. 
The  number  of  samples  were  found  to  be  sufficient  because  further  increase  in 
sample  size  did  not  significantly  effect  the  output  variable  distribution  and 
probability  limits.  [Ref.  45  :p.  8] 

Figure  60  shows  a  scattergram  of  the  state  vectors  which  define  the  entry  flight  path 

angle  [Ref.  45  :p.  7].   The  trajectories  were  then  simulated  from  the  entry  points  down 

to  the  breakup  altitude  for  each  randomly  selected  state  vector.  Figures  61  and  62  show 

the  results  of  a  nominal  run  used  to  determine  the  flight  characteristics  of  Skylab  [Ref. 

45: pp.  8-9].   From  Figure  62,  the  following  parameters  vs  altitude  were  presented: 


221 


Table  XVI:  ENTRY  STATE  VECTOR  /  ERROR  COVARIANCE  MATRIX 
[Ret.  45] 


Nominal  Entry  Point 

x  »  17,643,032  ft       V  =  -3238.58  ft/sec 

A 

y  =  -5,174,810  ft       Vy  =  21011.96  ft/sec 
z  =  -10,760,572  ft       Vz  -  -14693.59  ft/sec 


ECI  Error  Covariance  Matrix 


V      V      .V, 
x      y      z 


x  y      z 

x  5.76546E8 

y  -3.95374E9  2.72016E10 

z  2.79469E9  -1.92075E10  1.35675E10 

V  4.86448E6  -3.34635E7  2.36298E7  4.11679E4 

A 

V  -1.35955E6  9.34064E6  -6.59878E6  -1.14913E4  3.20968E3 

y 

V  -3.03260E6  2.08426E7  -1.47220E7  -2.56422E4  7.16010E3  1.59759E4 
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Figure  60:  Relative  Entry  Flight  Path  Angle  Scattergram 
[Ref.  45] 
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Figure  61:  Skylab  Nominal  Flight  Characteristics 
[Ref.  45] 
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Figure  62:  Downrange  Impact  Point  Variations 
[Ref.  45] 
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1 .  Dynamic  pressure 

2.  Mach  number 

3.  Load  factor 

4.  Velocity 

5.  Relative  flight  path  angle 

In  Figure  62,  the  275,000  ft  (84.8  km)  breakup  altitude  of  the  simulation  is  shown  along 
with  the  downrange  variations  for  several  debris  sizes  (denoted  by  ballistic  parameter  or 
coefficient).  The  lat/long  histogram  and  cumulative  probability  distribution  diagrams 
respectively  shown  in  Figures  63  and  64  [Ref.  45:pp.  11-12]  indicate  the  two  distinct 
regions  corresponding  to  the  largest  and  smallest  ballistic  coefficient  debris. 
Additionally,  these  diagrams  determine  the  lat/long  bounds  where  all  debris  is  likely  to 
fall.  [Ref.  45:pp.  7-12] 

The  downrange  and  crossrange  debris  impact  dispersion  distances  were 
calculated  using  the  lat/long  and  azimuth  angles  at  the  nominal  impact  point  by  the 
matrix  equation  [Ref.  45  :p.  12] 
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[Ref.  45] 
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Figure  64:  Longitude  Histogram  /  Cumulative  Probability  Distribution 
[Ref.  45] 
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N  =  nominal 

P  =  perturbed 

Figure  65  shows  the  downrange  and  crossrange  dispersion  distances  from  the  nominal 

impact  point  calculated  from  equation  (216)  [Ref.  45:p.l4].  As  shown  in  the  figure,  the 

crossrange  dispersions  were  relatively  small  as  compared  to  the  downrange  dispersions. 

In  order  to  determine  the  effect  of  the  breakup  altitude  on  the  dispersion  area, 

Monte  Carlo  simulations  were  conducted  at  various  altitudes  from  200,000  to  350,000 

ft  (61.7-107.9  km).    In  this  altitude  range,  the  reentry  object  experiences  large  thermal 

and  structural  loads.    Figure  66  shows  three-sigma  downrange  and  uprange  dispersion 

for  the  largest  and  smallest  debris  [Ref.  45 :p.  15]. 

This  figure  is  useful  for  estimating  the  total  footprint  dispersion  resulting  from 
breakup  of  smaller  pieces  at  higher  altitudes  and  heavier  pieces  at  lower  altitudes, 
thereby  accounting  for  the  uncertainty  in  breakup  altitude.  [Ref.  45:p.l4] 

In   their  conclusions,    the   authors   note   significant   footprint  dispersion 

variations  for  breakup  in  the  200,000  to  350,000  ft  altitude  region.   Furthermore,  they 

state  that  the  comprehensive  Monte  Carlo  analysis  is  very  useful  and  appropriate  in 

determining  impact  dispersion  areas  of  a  spacecraft  or  discarded  portions  of  a  launch 

vehicle.  [Ref.  45:p.  15] 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.      CONCLUSIONS 

The  primary  goal  of  this  thesis  was  to  identify  the  "state-of-the-art"  of  orbit-decay- 
induced  uncontrolled  reentry  and  impact  prediction  methods,  with  an  emphasis  on  the 
physics  of  the  final  few  revolutions  to  impact.  This  was  accomplished  through  a 
comprehensive  literature  survey  from  the  1950's  to  the  present  of  unclassified  military 
and  civil  databases.  The  survey  indicated  that  there  is  some  significant  foreign  work 
being  done  and  much  of  it  was  not  available  to  the  authors,  in  English  translation,  and 
thus  was  not  included  in  this  survey.  Also,  the  literature  survey  reflects  the  changing 
scientific  terminology  over  the  course  of  several  decades  and  it  is  especially  noticeable 
in  the  different  forms  that  the  common  variables  take  in  the  numerous  equations 
presented.  The  authors  did  not  make  any  attempt  to  use  the  standard  AIAA 
astrodynamics  nomenclature  or  to  standardize  the  equations  in  any  other  way. 

The  principal  conclusion  of  this  thesis  is  that  the  current  uncontrolled  reentry  and 
impact  prediction  methodology,  used  in  the  U.S.  and  abroad,  is  based  on  analysis  which 
is  30  or  more  years  old.  This  conclusion  is  based  on  the  fact  that  the  U.S.  method  takes 
its  roots  in  the  works  of  Brouwer  (1959)  and  Allen  and  Eggers  (1957),  and  that  the  U.S. 
method  is  the  accepted  international  standard,  as  shown  by  the  literature  survey. 

While  conducting  the  literature  search  dating  back  to  the  1950's,  the  authors 
noticed  a  definite  trend,  through  the  years,  in  the  focus  of  published  material  related  to 
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reentry.   During  the  timeframe  from  the  late  1950's  until  the  mid  1960's,  the  emphasis 

in  the  literature  was  primarily  devoted  to  understanding  and  describing  the  physics  of 

reentry,  specifically  as  it  pertains  to  controlled  ballistic  reentry  of  Mercury  and  missile 

type  vehicles.    With  the  exception  of  Sputnik  IV,  very  little  of  the  literature  surveyed 

during  this  timeframe  was  strictly  devoted  to  uncontrolled  reentry.  However,  during  the 

derivations  of  their  analytical  reentry  theories,  some  early  pioneers  such  as  Chapman  and 

Loh  investigated  the  shallow  orbit-decay-induced  reentry.    With  the  physics  of  reentry 

fairly  well  understood,  starting  in  the  late  1960's  and  continuing  into  the  1980's,  the 

emphasis  in  the  literature  had  shifted  to  controlled,  gliding  reentry  in  order  to  support 

the  launch  of  the  Space  Shuttle.     In  1965,   Eggers  and  Cohen  stated: 

Significant  advances  have  been  accomplished  in  the  science  and  technology 
appropriate  to  atmosphere  entry  of  spacecraft  during  the  eight  years  since  the 
launching  of  Sputnik  I.  This  progress  is  illustrated  by  the  successful  entry  from 
Earth  orbit  of  the  manned  Vostok,  Mercury,  Voskhod,  and  Gemini  spacecraft,  for 
which  the  problems  in  orbital  entry,  such  as  high  convective  heating  rate  and  load 
and  communication  blackout,  were  successfully  over  come... Although  the  first 
Apollo  vehicle  entry  has  yet  to  be  demonstrated... Much  engineering  work  for  this 
vehicle  remains  to  be  done;  however,  the  fundamental  research  activity  associated 
with  Apollo  is  being  reduced  in  favor  of  that  associated  with  missions  and  vehicles 
of  the  more  distant  future.  [Ref.  75:pp.  339-340] 

Uncontrolled  reentry  briefly  came  back  into  focus  in  the  late  1970's  and  early  1980's 

with  the  reentry  of  Cosmos-954  and  Skylab  and  continued  throughout  the  1980's  in 

varying  degrees  where  the  reentries  of  Cosmos- 1402,  1601  and  Salyut-7/Cosmos-1686 

continued  to  spark  a  flurry  of  literature  from  the  European  Space  Agency. 

The  literature  survey  of  recent  publications  shows  a  strong  reliance  on  work  done 

in  the  1950's  and  1960's  as  a  basis  for  "extensions"  or  "modifications"  of  pre-existing 
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methods  used  for  prediction  of  reentry  and  impact.  Although  reliance  on  pre-existing 
methods  is  not  in  and  of  itself  flawed,  the  modern  works  often  subscribe  to  the  same 
assumptions  which  the  original  authors  were  forced  to  make  because  of  hardware  or  data 
limitations  of  their  time.  This  serves  to  perpetuate  the  inherent  limitations  of  the 
method's  ranges  of  applicability  and  validity.  For  example,  it  has  been  shown  [Refs. 
14,16-17]  that  some  very  important  work  done  in  the  1950's,  such  as  Allen  and  Eggers 
[Ref.  65],  which  is  still  routinely  referenced  in  recent  publications,  is  inherently  flawed 
due  to  assumptions  made  under  conditions  of  little  or  no  data.  This  particular  reference 
is  especially  important  since  it  is  the  basis  for  most  of  the  modern  reentry  heating  work. 
It  becomes  even  more  significant  since  the  literature  survey  has  shown  the  coupled 
dependence  of  reentry  breakup  to  reentry  heating  and  dynamic  load  effects.  Reference 
[42]  describes  the  limitations  imposed  by  these  assumptions  through  mathematical  proofs. 
This  is  further  substantiated  by  reentry  breakup  observations  [Ref.  86]  which  indicate  that 
the  classical  convective  heat  transfer  equations,  when  applied  to  breakup  analysis, 
consistently  underestimate  reentry  survivability. 

Of  the  various  "extensions"  to  the  current  reentry  theory,  of  which  the  NORAD 
method  is  recognized  as  the  international  standard,  there  does  not  appear  to  be  any  one 
method  which  is  singularly  superior  to  the  others.  However,  it  is  the  opinion  of  the 
authors  that  the  ESA  FOCUS  program  merits  special  attention  for  further  research.  This 
conclusion  is  based  on  the  fact  that  this  program  contains  a  very  sophisticated  method  for 
dealing  with  the  drag  coefficient,  CD,  and  could  easily  incorporate  other  higher-order 
effects  occurring  in  reentry,  such  as  attitude  generated  lift,  rotation  of  the  orbital  plane 
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(due  to  atmosphere  rotation  at  attitudes  <  60  km),  linear  variations  in  the  ballistic 
coefficient  in  the  final  few  days  of  orbit  and  others  as  mentioned  previously  in  this  thesis. 

Numerous  works  cited  in  the  literature  survey  also  describe  computing  power  and 
computer  time  constraints  as  critical  parameters  of  the  problem.  These  hardware 
limitations  then  forced  one  or  more  simplifying  assumptions  to  be  made  in  order  to  deal 
with  the  technology  limitations  of  the  day.  Subsequently,  much  of  the  literature  survey 
shows  work  which  has  limited  applicability  for  the  problem  of  concern. 

Additionally,  there  is  a  general  lack  of  consensus  in  the  literature  as  it  pertains  to 
orbit-decay-induced  uncontrolled  reentry.  The  very  definition  of  "reentry"  remains  vague 
and  often  is  defined  differently  by  investigators  attempting  to  do  similar  research.  Many 
investigators  have  examined  various  aspects  of  the  problem;  however,  when  surveying 
the  literature,  a  common  approach  and  standardized  starting  point  for  solving  this  very 
dynamic  problem  is  not  apparent.  A  good  example  to  illustrate  this  point  is  the  effect 
of  angle  of  attack,  a,  on  the  ability  to  accurately  predict  uncontrolled  reentry  and  impact. 
The  ability  to  properly  model  the  uncertainties  in  the  reentry  body  configuration  such  as 
changes  in  area,  mass,  and  attitude  is  a  very  difficult  problem.  As  a  result  of  these 
changes,  three-dimensional  angle  of  attack  becomes  an  exceedingly  difficult  parameter 
to  characterize  especially  when  the  reentry  body  is  undergoing  rapid  configuration 
changes  due  to  ablation  and  structural  deformation.  However,  in  order  to  properly  model 
the  aerodynamic  coefficients  for  lift  and  drag,  angle  of  attack  must  be  considered. 
Specifically,  lift  as  a  function  of  angle  of  attack  not  only  affects  the  trajectory  (location 
of  impact),  but  it  also  affects  the  altitude  of  breakup  (survivability  and  dispersion  area) 
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since  it  reduces  the  thermal  and  structural  loads  on  the  vehicle.  However,  in  most  cases 
where  angle  of  attack  is  investigated,  this  coupled  effect  is  not  mentioned  nor  pursued. 
Even  more  likely  in  the  literature,  the  reentry  body  is  modeled  as  a  point  mass  where 
angle  of  attack  is  not  considered  or  lift  is  assumed  to  be  zero  or  negligible. 

Another  observation  from  the  survey  is  that  no  comprehensive  sensitivity  and  error 
analysis  has  been  conducted  in  order  to  determine  quantitatively  the  effects  of  critical 
parameters/variables  on  impact  prediction.  For  example,  numerous  sensitivity  and  error 
analysis  studies  have  been  conducted  on  atmospheric  density  models  uncoupled  from 
other  critical  parameters.  However,  a  coupled  analysis  of  all  of  the  critical  parameters 
could  determine  which  variables  contribute  most  significantly  to  the  overall  accuracy  of 
the  various  impact  prediction  methods. 

In  addition  to  the  author's  observations  and  conclusions  described  in  the  previous 
paragraphs,  a  summary  of  specific  conclusions  from  the  literature  survey  are: 


1.  There  is  no  clear  definition  of  when  reentry  occurs.  Typically  120  km  or  an 
orbital  period  of  87.5  minutes  or  less  is  used  in  the  various  models  to  define  the 
start  of  reentry.  However,  observations  and  post  flight  analysis  indicate  a  range 
of  altitudes  where  "reentry"  actually  occurs. 

2.  There  is  a  lack  of  observation  data  in  the  reentry  regime  due  to  a  lack  of  global 
sensor  coverage.  This  lack  of  data  significantly  adds  to  the  uncertainty  of  the 
problem. 

3.  The  current  deterministic  dynamics  model  appears  to  be  inadequate  for 
processing  the  true  physics  of  reentry. 

4.  The  lack  of  observational  data  coupled  with  the  inadequate  knowledge  of  the 
true  physics  of  uncontrolled  reentry  significantly  increases  the  uncertainty  of  the 
problem. 
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5.  Reentry  is  most  likely  to  occur  within  several  tenths  of  a  degree  of  the  minimum 
altitude  (height  above  the  ground)  point  in  the  orbit,  which  is  not  necessarily  the 
point  of  perigee.  [Ref.  99] 

6.  The  uncontrolled  reentry  is  three  times  more  likely  to  occur  in  a  latitude  band 
equaling  the  inclination  than  in  an  equatorial  band  (0  to  10  degrees),  with  the 
exception  of  polar  orbits  where  an  equatorial  reentry  is  most  likely  to  occur. 
[Ref.  99] 

7.  The  impact  location  (downrange  and  crosstrack)  is  affected  by:  [Refs.  45,101] 

(a)  the  rotation  of  the  atmosphere,  starting  at  an  altitude  of  30-60  km 

(b)  the  debris  ballistic  coefficient 

(c)  the  entry  dispersions  (velocity  and  position  errors  at  the  entry  altitude) 

(d)  the  relative  flight  path  angle 

(e)  the  breakup  altitude 

(f)  the   environmental   conditions   such   as   wind   direction/magnitude  and 
atmospheric  density 

8.  There  appears  to  be  linear  variations  in  the  ballistic  coefficient  in  the  final  few 
days  prior  to  reentry.  [Ref.  92] 


B.      RECOMMENDATIONS 

Based  on  the  above  conclusions,  the  following  recommendations  are  proposed: 

1.  As  mentioned  in  Chapter  I,  there  was  no  attempt  made  to  standardize  the 
nomenclature  and  variables  of  the  equations  presented  in  this  thesis,  in 
accordance  with  AIAA  astrodynamics  standards  or  any  others.  Based  on  the 
literature  survey,  there  clearly  exists  a  need  for  such  standards,  especially  when 
dealing  with  work  spanning  the  course  of  several  decades. 
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2.  This  literature  survey  was  conducted  through  unclassified  sources,  written  in 
English,  only.  A  survey  of  classified  and  foreign  language  sources  should  be 
conducted. 

3.  The  primary  focus  of  this  thesis  was  a  thorough  survey  of  the  physics  of 
uncontrolled  reentry.  Throughout  the  survey,  the  authors  noted  several 
operational  issues  which  should  be  addressed  in  follow-on  research.  These  issues 
include  human  factors  (orbit  analyst  experience),  sensor  bias  and  coverage,  and 
improved  international  information  exchanges. 

4.  Several  reentry /impact  prediction  methods  were  presented  (LIFETIME,  FOCUS, 
and  NORAD)  in  this  thesis.  A  side  by  side  comparison  of  these  methods  using 
historical  data  should  be  conducted  in  order  to  compare  and  contrast  their 
overall  accuracies  and  efficiencies. 

5.  A  sensitivity  and  error  analysis  of  critical  parameters  including  atmospheric 
density,  aerodynamic  coefficients,  initial  conditions  (vehicle  area,  mass,  attitude, 
position  and  velocity),  thermal  and  dynamic  loads,  and  breakup  phenomena 
should  be  conducted  in  order  to  improve  the  current  dynamics  model  and  to 
focus  future  research  efforts. 

6.  In  view  of  the  lack  of  global  sensor  coverage,  especially  in  the  southern 
hemisphere,  high  priority  reentry  events  may  be  better  predicted  through  a  joint 
cooperative  effort  using  Navy  and  Air  Force  assets  as  mobile  supplemental 
observation  and  tracking  stations. 

7.  Since  uncontrolled  reentry  is  characterized  by  numerous  rapidly  changing 
critical  parameters,  which  are  ill-defined,  stochastic  and  statistical  methods 
should  be  applied  to  current  reentry  models  to  better  analyze  the  sensitivity  of 
the  various  uncertainties  associated  with  this  problem.  This  could  serve  as  an 
"operational  tool"  to  help  improve  the  prediction  accuracies  of  the  current 
models  while  the  physics  of  the  reentry  process  is  being  more  thoroughly 
investigated. 

8.  The  authors  strongly  recommend  future  cooperative  research  between  the  Naval 
Postgraduate  School  and  the  AFSPACECOM/USSPACECOM  in  order  to  solve 
this  largely  unknown  and  highly  dynamic  process. 
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